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-1- 

Toepassing van a priori informatie in de "stacking" procedure, gebruikt by de analyse van 
oppervlaktegolven, leidt tot verbetering van de identificatie van individuele boventonen van 
opperviaktegolven, niet tot een nauwkeuriger bepaling van fase- of groepssnelheid van goed 
geidentificeerde boventonen. 


-2- 
Metingen van fasesnelheden van grond- en boventonen van Love golven in West-Europa 
zijn tot dusver onbetrouwbaar. 
Cara, M., Nercessian, A. and Nolet, G., 1980. 
Geophys. J. R. Astron. Soc., 61 , 459-478. 
Dit proefschrift 


eels 
Een mogelijke inversie van de dichtheid op een diepte van ca. 220 km (Nolet, 1977) wordt 
niet gevonden onder het Westeuropees platform. 

Nolet, G., 1977. 

Z. Geophys., 43, 265-286. 

Cara et al., 1980 

Dit proefschrift 


4. 
Het ontbreken van een lage snelheidslaag voor "S"- golven en de aanwijzingen voor de 
aanwezigheid van een lage dichtheidslaag op een diepte tussen 200 en 350 km diepte onder 
het Scandinavische schild zijn in overeenstemming met de tectosphere hypothese. 

Jordan, T. H., 1978. 

Nature., 274, 544-548. 

Dit proefschrift 


-5- 

By de interpretatie van aardmodellen, bestaande uit gegevens over de "S"- snelheid en 
dichtheid, in termen van mineralogie verdient het de voorkeur ontbrekende informatie 
omtrent elastische parameters van mantel-mineralen als onbekende te beschouwen in plaats 
van te vertrouwen op extrapolatie van bekende gegevens. 


7 Be 
Micro-electronica wordt sneller ontwikkeld dan getest. 


Es 
De hoeveelheid problemen die onstaat bij niet-gestandaardiseerde douanehandelingen is 
evenredig aan het aantal verdragen tussen de landen in kwestie. 


-8- 
Het aantal projecten dat gestimuleerd wordt door de Europese Gemeenschap is geen maat 
voor de de bereidheid tot samenwerking van de diverse landen. 


-9- 

Met het verdwijnen van studierichtingen door toedoen van de huidige bezuinigingsgolf rijst 
de vraag of de aanduiding Universiteit in de betekenis van "het geheel der wetenschappen" 
nog wel op zijn plaats is. 


-10- 
De kritiek dat popmuziek inhoudsloos zou zijn wordt weerlegd door teksten als: "brown 
shoes don’t make it" (F. Zappa). 


Utrecht, 13 mei 1987 B. Dost 


Stellingen behorende bij het proefschrift: “The NARS array”. 
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SUMMARY 


Due to the rapid development of portable, digital seismographs it has recently become pos- 
sible in global seismology to install and operate a large scale temporary array of seismic sta- 
tions. This thesis describes the design and operation of the first experiment of this kind: the 
Network of Autonomous Recording Stations (NARS). The aim of NARS is to measure the 
dispersion of surface waves, especially higher modes, and to obtain through linearized 
inversion a shear-velocity and density model for the upper mantle under the West European 
platform. The existence of such a regional model opens the possibility to test different 
current hypotheses about the structure and composition of the upper mantle. Inversion of 
previous measurements of the dispersion of surface waves in Western Europe resulted in an 
average model for the West European platform including the Scandinavian shield. Combi- 
nation of the measurements presented in this thesis with the older dataset opens the possibil- 
ity to investigate differences between the shield and platform structure. 


In the first chapter the design and a technical description of a NARS station is given. Furth- 
ermore NARS installation and operation is discussed. Inadequate sampling due to instru- 
ment failures and absence of strong Japanese/Kurile events at intermediate depth, which 
excite most efficiently higher mode surface waves, are the main sources of error in the 
presented dispersion measurements. These facts motivated us to develop methods to 
enhance the signal to noise ratio in the stacking technique applied to separate higher modes 
of surface waves. 

An inventory of possibile methods to obtain this enhancement followed by a test of the pro- 
posed methods on a synthetic dataset is the subject of chapter 2. Introduction of a priori 
information in the stacking procedure increases the signal to noise ratio and facilitates the 
identification of a mode, but does not improve the average error in the phase velocity deter- 
mination of higher mode Rayleigh waves. However, the use of the CLEAN-algorithm in 
removing disturbing aliasing lobes from well excited modes, is demonstrated to produce a 
bias in phase velocity measurements at frequencies greater than 35 mHz. 

Chapter 3 describes the dataset that has been collected by the NARS array and could be 
used in the stacking procedure. In addition phase velocity measurements for fundamental 
and higher Rayleigh and Love waves are presented. Higher mode Rayleigh wave phase 
velocities could be obtained in the frequencyband of 20-70 mHz, an extension of 30 mHz to 
the high frequency part with respect to earlier measurements, and fundamental mode phase 
velocities from 6-50 mHz. Because of the very limited amount of reliable transverse com- 
ponents, the Love wave phase velocity dataset has been rejected for further use in the inver- 
sion procedure. An independent dataset, consisting of Sn traveltimes from selected Euro- 
pean events and recorded by European stations, is presented in order to constrain the shear 
velocity in the uppermost mantle. 


Inversion of the phase velocity dataset is the subject of chapter 4. A detailed shear velocity 
and density model for the West European platform has been obtained. Major features are a 
high velocity lid from 80-140 km depth, followed by a pronounced low velocity zone 
between 160 and 220 km depth. The transition zone (400-650 km depth) is characterized by 
high shear velocity values with respect to a reference Earth model (PREM). Density shows 
a positive gradient coinciding with the high velocity lid, followed by a constant value down 
to the 400 km discontinuity. No low density zone is resolved. Combination of Nolet’s 
(1977) higher mode phase velocity dataset for the West European platform and Scandina- 
vian shield and the present phase velocities for the platform structure resulted in a model for 
the Scandinavian shield. Remarkable features are the absence of a low velocity layer, low 
velocities compared to the platform model in the upper part of the transition zone and low 
densities between 200 and 350 km depth. 

The final chapter 5 gives an interpretation of the models discussed in chapter 4 in terms of 
mantle mineralogy. Interpretation of the platform model in terms of pyrolite or piclogite is 
not conclusive. Density, especially around the 400 km discontinuity, indicates that the 
pyrolite model is to be preferred. The high velocity layer associated with an increase in 
density is interpreted as an eclogite layer. The existence of this layer can be explained by 
the lithospheric doubling concept. Difference in structure between platform and shield, pos- 
sibly down to 500 km depth, strongly supports the tectosphere hypothesis. 


SAMENVATTING 


Dankzij de snelle ontwikkeling van draagbare, digitale seismografen is in de globale 
seismologie de mogelijkheid geschapen om een tijdelijk netwerk van seismische stations te 
installeren en daarmee een experiment op grote schaal uit te voeren. 

Dit proefschrift beschrijft de opzet en uitvoering van het eerste experiment dat met een der- 
gelijk tijdelijk netwerk is uitgevoerd: het Netwerk van Autonoom Registrerende Stations 
(NARS). 

Het doel van het hier beschreven experiment is om een gedetailleerd S-snelheids en 
dichtheids model voor het Westeuropees platform te verkrijgen door middel van inversie 
van dispersiemetingen van grondtoon en (met name) boventonen van oppervlaktegolven. 
Een zo verkregen regionaal model voor de bovenmantel kan gebruikt worden om heersende 
theorieen omtrent structuur en compositie van de bovenmantel te toetsen. Inversie van eer- 
dere metingen van de dispersie van oppervlaktegolven in West-Europa leverde een gemid- 
deld model op van het Westeuropees platform en het Scandinavisch schild. Combinatie van 
deze eerdere metingen met de huidige geeft de mogelijkheid om verschillen tussen beide 
tegio’s te onderzoeken. 


De opzet, het technisch ontwerp en de uitvoering van NARS is beschreven in hoofdstuk 1. 
Het ontbreken van sterke aardbevingen op intermediaire diepte, die een optimale excitatie 
van de boventonen geven, en het frequent ontbreken van data uit NARS-stations die tijdelijk 
defect geraakt waren, zijn de belangrijkste oorzaken van fouten in de dispersie metingen. 
Dit heeft geleid tot de ontwikkeling van methoden om de signaal-ruis verhouding te vergro- 
ten in het stacking-algorithme dat gebruikt wordt om boventonen van elkaar te scheiden. 
Een inventarisatie van mogelijkheden hiertoe en een test op een synthetische dataset is het 
onderwerp van hoofdstuk 2. Introduktie van a priori kennis in de analyse geeft een 
belangrijke verhoging van de signaal-ruis verhouding, waardoor verschillende boventonen 
gemakkelijker te identificeren zijn, maar verkleint niet de gemiddelde fout in de metingen. 
Aangetoond wordt dat het CLEAN-algorithme, dat iteratief het effect van een niet optimale 
Sstationsbedekking corrigeert, fouten introduceert voor frequenties groter dan 35 mHz. 

Een overzicht van de door NARS geregistreerde aardbevingen die gebruikt konden worden 
in de analyse staat in hoofdstuk 3. Ditzelfde hoofdstuk behandelt tevens de fasesnelheids 
metingen van grond- en boventonen van Rayleigh en Love golven. Fasesnelheden van Ray- 
leighgolven worden gepresenteerd in de frequentieband van 20-70 mHz voor de boventonen 
en van 6-50 mHz voor de grondtoon. Dit is een extensie van 30 mHz naar hogere frequen- 
ties met betrekking tot eerdere metingen in West-Europa. Vanwege de geringe hoeveelheid 
betrouwbare transversale componenten en de grote spreiding in de gemeten fasesnelheden 
van Love-golven, zijn deze niet meegenomen in de inversie procedure. Een onafhankelijke 
dataset bestaande uit looptijden van Sn-golven, afkomstig van Europese aardbevingen en 


Ii 


geregistreerd door Europese stations, is geintroduceerd om de S-snelheid in het bovenste 
deel van de bovenmante! vast te leggen. 

Het onderwerp van hoofdstuk 4 is de inversie van de in hoofdstuk 3 gepresenteerde dataset. 
Een gemiddeld model voor het Westeuropees platform wordt gekarakteriseerd door een 
zone met hoge snelheden tussen 80 en 140 km diepte, gevolgd door een geprononceerde 
lage snelheidslaag tussen 160 en 220 km diepte. De overgangszone (400-650 km diepte) 
bevat hoge snelheden ten opzichte van een referentie aardmodel (PREM). De dichtheid 
geeft een positieve gradient te zien die samenvalt met de zone van hoge snelheid, gevolgd 
door een gebied met een nagenoeg konstante dichtheid tot 400 km diepte. Er is geen 
aanwijzing voor een lage dichtheidslaag. Vergelijking met data van Nolet (1977) voor plat- 
form en schild tesamen resulteerde in een model voor het Scandinavisch schild. Opvallend 
is het ontbreken van een lage snelheidslaag, lagere snelheden ten opzichte van het platform 
model in het bovenste deel van de transitiezone en lagere waarden in dichtheid van ca. 200 
tot 350 km diepte. 

In hoofdstuk 5 wordt een mineralogische interpretatie van deze snelheids en dichtheidsmo- 
dellen gegeven. Er kan geen uitsluitsel gegeven worden omtrent de vraag of het platform 
model een pyrolitische ofwel een piclogitische samenstelling heeft. De dichtheid, en wel 
speciaal de dichtheidssprong bij de 400 km discontinuiteit, geeft echter aanwijzingen dat het 
pyroliet model de voorkeur geniet. De hoge snelheidszone van 80 tot 140 km diepte, die 
gepaard gaat met een verhoging van de dichtheid, wordt geinterpreteerd als een eclogiet- 
laag. De aanwezigheid van deze laag kan verklaard worden door middel van de lithosfeer- 
verdubbelings hypothese. Verschillen in structuur tussen platform en schild, mogelijk tot op 
500 km diepte, zijn in overeenstemming met de tectosfeer hypothese. 


IV 
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INTRODUCTION 


Layered models of the deep upper mantle structure summarize information on the 
average structure of the sampled region. Although there is a growing tendency in seismol- 
ogy to shift focus to laterally heterogeneous models and their implications for the dynamic 
and kinematic processes of plate tectonics, simplified models with horizontal layers shall 
continue to be useful for several reasons. First of all they may provide much-needed evi- 
dence for petrologic models of the upper mantle, especially in relatively undisturbed areas 
like the West European platform. Secondly, the density in the upper mantle is still badly 
constrained, and even the simple layered models coming from surface wave analyses may 
constitute an important step forward. Thirdly these layered models are needed as a starting 
model in more detailed studies, e.g. in linearized or non-linear waveform fitting (Lerner- 
Lam & Jordan, 1983; Nolet et al., 1986d), or as a description of the background medium in 
body wave (Malin & Phinney, 1985) and surface wave scattering calculations (Snieder, 
1986). For the West European platform no deep layered model is available at present. A 
study by Mayer-Rosa and Mueller (1973) provided such a model for south West Europe 
derived from fundamental mode Rayleigh wave phase velocities and P and S delay times. 
The deep mantle models that have been derived from inversion of fundamental and higher 
mode Rayleigh wave phase velocities (Nolet, 1977) and Love and Rayleigh wave funda- 
mental and higher mode phase velocities (Cara et al, 1980) averaged over platform and 
shield. Recently tomographic studies using longitudinal body waves (Spakman, 1986) give 
detailed pictures of the upper mantle structure. Shear wave tomography, however, has not 
yet produced accurate results. 

Surface wave dispersion is mainly dependent on the shear wave velocity and density 
structure of the traversed medium, if one assumes the medium to be isotropic. Nolet (1977) 
showed that linearized inversion of higher mode phase velocities, taking averaged measure- 
ments from 7 events, for shear velocity alone did not produce accurate earth models. Addi- 
tion of the density in the inversion resulted in the appearence of a low density layer starting 
at 190 km depth. However, resolution was not sufficient to constrain this low density layer 
and its existence was questioned by Cara et al. (1980) who took a subset of Nolet’s data and 
added a Love wave dataset. The main reasons for this bad resolution are: 


[a] In order to obtain data from a large enough number of stations, both Nolet (1977) and 
Cara et al. (1980) average over a heterogeneous structure; 


[b] both used analogue WWSSN recordings are used with digitization errors increasing 
with frequency; 


[c] the large separation between WWSSN stations leads to disturbing spatial aliasing 
effects so that phase velocity measurements could only be obtained at frequencies < 
40 mHz. 


If one wishes to measure higher-mode surface wave dispersion without the incor- 
poration of a priori knowledge in the form of a starting model, stacking techniques should 
be applied. The higher mode stacking techniques developed by Nolet (1975,1977) and Cara 
(1978) require the existence of a dense, quasi-linear array of seismic stations. In an analysis 
of measurement errors in higher-mode phase velocity determinations, Nolet and Panza 
(1976) find that the largest source of errors is the interference of different modes with 
roughly the same group velocity combined with an inadequate array-response. Precision 
largely depends on the array span, which determines the width of the mainlobe of the array 
response function, and station density. Addition of temporary digital stations to the existing 
WWSSN network should improve the higher-mode dataset considerably. 

Based on these conclusions and the wish to sample a geological unit like the West 
European platform, the Network of Autonomous Recording Stations (NARS; Nolet and 
Vlaar, 1982) was initiated. This thesis describes the design and set-up of the network, 
which consists of 14 digital, portable, 3 component stations located in a linear array from 
south Sweden to southern Spain. This gives a total array length of 2600 km, or an average 
Station separation of 200 km. The array acts as an antenna directed towards the 
Kuriles/Japan region and was designed to deliver high quality data and a higher mode phase 
velocity dataset in the frequency range of 20-70 mHz. 

The final aim of the research described in this thesis is to obtain a detailed picture of 
the average West European platform structure by measuring the dispersion of higher mode 
surface waves and increase the vertical resolution of shear velocity and density with respect 
to earlier studies. This enables us to test different hypotheses regarding the constitution of 
the upper mantle. 


NARS is the Dutch contribution to the European Geotraverse (EGT). 


PART 1: THE NETWORK 
1. DESCRIPTION OF THE NARS ARRAY 


For many years seismological observatories have been equipped with non-mobile 
seismometers and analogue recording. Maintenance could only be carried out by qualified 
personnel. A common seismic station consisted of a seismometer- galvanometer system, 
producing an output on photographic paper or ink-recorders. With the development of port- 
able seismometers and, more recently, the replacement of galvanometers by an electronic 
amplifier-filter system (Wielandt and Mitronovas, 1976) and the introduction of commer- 
cially produced digital event recorders with triggered recording, it became possible to carry 
out a large scale seismic experiment using portable stations. In this chapter we will describe 
the NARS system design and give detailed technical information about the instrumentation. 
Furthermore we will discuss the operation of the array and the data processing in Utrecht. 
Because of recent developments in developing a large broad-band portable network in the 
US CRIS, 1984b) and in Europe (Nolet et al., 1986c) a brief overview of our operational 
experience with the NARS instrumentation is included. For the upgrading of the WWSSN 
network (IRIS, 1984a) operational experience from the IDA network has recently been 
reviewed (Agnew et al., 1986). Parts of this chapter have been published in Dost et al. 
(1984) and Nolet et al. (1986a). 


1.1. NARS SYSTEM DESIGN 


The NARS instrumentation was selected with the following requirements in mind: 


{a] The array should record large earthquakes with a good signal to noise ratio at 3 com- 
ponents. Since the aim of the network is to study Earth structure, not seismicity, we 
are not interested in events of small magnitude (mm, <6); 


{b] The frequency band should include surface waves in the period range between 10 and 
200 mHz, and give broad-band ( < 2 Hz) recordings of body waves; 


{c] Digital recording and a large dynamic range is necessary to allow for sophisticated 
data analysis techniques. Seismic signals of interest have a dynamic range of at least 
120 dB; 


{d] Mobility and ease in maintenance is required to make the network multi-purpose and 
low-cost. 


From [b] a lower limit for the Nyquist frequency was determined (f= 4 Hz) to avoid alias- 
ing. This means a minimum sampling rate of 8 samples per second. Continuous digital 


recording at this sampling rate, however, implies approximately 1.1 Mbyte of data per chan- 
nel per day. No event recorder was commercially available in 1981 that could store these 
large data volumes if we assume that point [d] requires station servicing at most once every 
week. In view of this we opted for triggered recording and a sampling rate of 125 msec. 
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Figure 1.1.1 NARS station configuration 


NARS station configuration is shown in figure 1.1.1. The central element of each station is 
formed by a Kinemetrics PDR-2 datalogger and 3 long period seismometers: one Teledyne 
Geotech SL 210 (vertical) and two SL-220 (horizontals). Matching the output of these 
seismometers to the PDR-2 necessitated the introduction of a pre-amplifier to take full 
advantage of the dynamic range of the PDR-2 (114 dB). The seismometers can be tuned to 
an eigenperiod of 10 to 30 seconds, but become increasingly unstable at longer periods. This 
made us decide to tune the seismometers to a 12 seconds eigenperiod and to introduce a 
response shaping filter after the pre-amplifier to enlarge the frequency band towards low 
frequencies. This pre-amplifier-filter was designed and built in the department of theoretical 
geophysics together with Arie van Wettum, who also built a time receiver for regular cali- 
bration of the internal PDR-2 quartz clock. 


1.2 NARS INSTRUMENTATION 
1.2.1. Seismometers 


The history of the development of portable long period seismometers shows a large 
effort to diminish the length of the spring of the vertical sensor and weight of the mass. 
With the introduction of the zero-initial length spring the portable seismometer could be 
realized. Both vertical and horizontal seismometers can be tuned to a specific eigenperiod 
by tilting the base of the instruments. This tilting results in a change of the gravity com- 
ponent in the plane of motion by cos 6 , where @ is the angle between the boom and the 
vertical (see e.g. Aki & Richards, 1980 pp 482-485). If the seismometer base is placed in a 


horizontal position (@=5) one can, theoretically, acquire an infinite eigenperiod. In prac- 


tice there is a finite range of eigenperiods where the instruments are stable. For NARS this 
region is between 10 and 30 seconds. We can write the relation between eigenperiod and 
equivalent mathematical spring length as: 


lovin 
T=2n (————_ 1. 
TS ease a) 


where / denotes spring length [m], g gravity constant [ms~*] and T the eigenperiod [s]. For 
10 seconds eigenperiod the equivalent spring length (//cos 6) is 25 m, but for 30 seconds the 
equivalent spring length is 225 m, Since the seismometer reacts to tilts of the base and 
extension of the spring, e.g. due to temperature fluctuations, as if the dimensions of the 
instrument are equal to the mathematical equivalent lengths (Willmore, 1979), one has to 
decide to tune the seismometers to a low eigenperiod for sake of stability, in our case 12.0 
seconds, 

The seismometer transferfunction (in Vm7') as function of complex circular frequency 


s=i @, defined as a response to a signal e*" is given by: 
3 3 
ss ss , 2\1/2 
T(s)= Se 21. 2 = Woh + 11-15) 
S?+2s oho?  (S-z;)(S-zq)’ 
(p= 0.524Hz hg = 1. (1.2) 


S is the main coil constant in Vm~!/Hz and is specified by the manufacturer to be 
$=90+4 Vm7'/Hz for the vertical instrument (SL-210), and $=93+4 Vm7'/Hz for the hor- 
izontal instrument (SL-220). 
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Figure 1.2.1 Pre-amplifier/filter circuit diagram. For explanation see text. 


The pre-amplifier-filter can be divided into three parts (see figure 1.2.1): 
[a] Pre-Amplifier, 

[b] Response shaping filter, 

{c] Additional amplification. 

In order to take full advantage of the large dynamic range of the PDR-2 (114 dB) a 
pre-amplifier was needed. The gain of this pre-amplifier is chosen in such a way that clip- 
ping occurs for an Italian event of magnitude 7 or a magnitude 8 teleseismic event. This 
gain factor has a value of approximately 40. Introduction of the shape response filter, how- 
ever, (see figure 1.2.3) decreases the signal level for frequencies higher than 0.1 Hz with a 


factor of about 65. A total gain of approximately 2600 should therefore be applied. If we 
apply this gain before filtering, the voltage after amplification can exceed the voltage of the 


power supply of the chopper amplifier (£15V ). This will introduce distortion of the signal. 
For this reason we restricted the amplifier gain to 510 before and 5 after filtering. A 
detailed description of the system, illustrated by figure 1.2.1, is as follows; 


[a] Pre-amplifier 


A central position is taken by the Analogic MP 221 chopper amplifier, with a specified noise 
level from DC to 1 Hz of less than 100 nV peak to peak and a drift less than 50nV/°C. Ry 
takes care of the critical damping of the seismometer. R, and R. determine the gain: 
R,+R2 
Ry 


gain = 


and R is used to diminish the DC offset. An offset error voltage is generated by the bias 
current flowing through the summing impedance of R, and R,. If we choose: 


R,R2 


a RR 


it compensates the offset, provided the source impedance is much less than R,. C, is a 
compensation capacity, used to limit the bandwidth of the system. This is important because 
output noise generated by the amplifier varies as the square root of the bandwidth. Its value 
has been calculated from the empirical formula: 


_ 30 100 
gain’ f 
where f is the cut-off frequency in Hz. This equation is determined by the manufacturer. 


Our choise of C,, see table 1.1, is based on a cut-off frequency of 15-20 Hz (in reality f= 
17.8 Hz). 


[WF] 


¢c 


[b] Response shaping filter 


In Appendix A the transfer function of a general second order filter is shown to be: 
Uy YiY2 


6 ee 1.3 
U; YiYsY,(V)+Y2+Y3) (1.3) 


U,, and U; denote respectively the output and input voltage and Y; admittance. If we iden- 
tify Y, and Y, with Rj! (see figure 1.2.1) and Y, and Y, with sC one obtains the well- 
known formula of a biquadratic low-pas filter: 


U, Wo 1 


U; S742 Moh tog ” ee RAC 


Incorporation of resistance Rs; in Y3 and Y, results in an additional high pass filter. 


sC 


In this case Y3 = Y, = ———— and: 
OR AERC 


U, 42s Whyte = o? 


U; ~—-s?42s@h;+m? § we 


ee ee ee 


where @, = RL’ 
5 


1 
(RatRs5)C’ 


(1.4) 


(1.5) 


The total response of the system consisting of seismometer and shaping filter is equal to the 
product of (1.2) and (1.5). If we choose Wp in (1.5) to be identical with the eigenfrequency 
of the seismometer ( Wp in (1.2)) we can change the eigenfrequency of the seismometer to 


any desired value @. 


{c] Additional amplification 


Instead of a chopper amplifier we use for the second amplifying stage an OPOS 
operational amplifier. The offset problem is handled in the same manner as in [a]. C, is 
introduced to cut off high frequencies, far outside the passband of the total system. Its value 


can be determined by: 


1 
a2 2n Ro f [FJ 


f is the cut-off frequency in Hz, resistance in Ohm. The choise of C,, displayed in table 1.1, 


implies a cut-off frequency of 39 Hz. 


Table 1.1: Component values 

R, 110kQ R> 56 kQ 
R3 110 kQ Rg, 10000 kQ 
Rs 1400 kQ Re 20 kQ 
R, 6.8 kQ Rg 5.1kQ 
Ro 27.2 kQ Rio 25 kQ 
Rg 3.8 kQ 

os 1400 nF Cz 330 nF 
C,; 150 nF 


1.2.3. Digital event recorder 


The PDR-2 digital event recorder is a low power multi-microcomputer based 

seismic data acquisition system. It operates with a triggering mechanism based on Short 
Term Average (STA) and Long Term Average (LTA) determination of the input signal. 
The parameter settings used are shown in table 1.2. 
A Kinemetrics CTU-1 play back unit is used for station servicing in the field and for a rough 
impression on paper of the data on cassette tapes that arrive in Utrecht. To inspect data in 
the field for e.g. parity errors a Tandy TRS-80 model 100 portable microcomputer is used. 
Data can be transmitted from the PDR-2 to the Tandy computer in the field at 9600 Baud. 
To accomplish this the commands from the TRS-80 to the PDR-2 are given through an RS- 
232 interface, but the data from the PDR-2 to the TRS-80 are received via the cassette 
recorder entrance. 


Table 1.2: PDR-2 parameter setting and general features 


3 channel recording 
125 msec sample period 
Trigger: 
STA/LTA ratio; no LTA update 
Only effective at vertical component 
STA = 45s; LTA = 128 s 
STA/LTA treshold = 4; Post event treshold = 2 
Post event time = 8*128 s 
All channels have auto gain ranging. 
Anti alias filter (2nd order Butterworth): -3dB point = 1.0 Hz 
Sample: 12 bit mantisse and 3 bit exponent 
Total number of samples per tape: 920.000 ( app. 10 hours of data) 
Recording : 
Standard 300 ft, 1667 bpi digital cassettes 
Blocks of 3072 samples (=6144 bytes) 
Header 114 bytes 


1.2.4. Time receiver 


Since our sampling rate is 125 msec, the clock precision needed is 63 msec or less. 
The PDR-2 can be equipped with a GOES-satellite time receiver, but signals from this satel- 
lite can not be received in Europe. We therefore had to build a time receiver ourselves and 
used a design published in Elektuur (Anonymus, 1981) that only needed a minor 
modification in order to synchronize the internal PDR-2 quartz clock. The role of this 


instrument in the total system is shown in figure 1.2.2, By means of an active antenna a 
coded time signal emitted by the DCF-77 station in Frankfurt is received. 
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Figure 1.2.2 Block diagram of the function of the time-receiver in connection 
with the PDR-2. 


At stations far away from Frankfurt an extra loop antenna is used, consisting of a coil 
(diameter of 107 cm, 78 windings) and capacitor (5-7 pF), tuned to 77.5 kHz. This 
increases the signal to noise ratio. The signal is decoded in a microprocessor unit and the 
time is shown on a display. With a modification of the PDR-2 we succeeded in establishing 
an automatical adjustment of the PDR-2 clock at pre-programmed time intervals. This pro- 
cedure synchronizes the clock to the nearest second; in consequence the clock drift must not 
be more than 0.5 seconds between subsequent synchronizations. To check for larger drift a 
minutemark is superposed on the seismic signal, on the East-West channel, during the first 
two minutes before each triggering. Since the first two recorded minutes of each event con- 
sist only of microseismic noise, no signal of interest is disturbed by this procedure. In prac- 
tice the drift of the internal PDR-2 clock can be more than 63 msec per week. Therefore the 
automatic time calibration is done twice a day (at 23.00 and 4.00 GMT) if reception is good 
and no earthquake is being recorded. 


1.2.5. System response 


The total system response is equal to the product of the individual responses of the 


seismometer, pre-amplifier/filter and anti-alias filter. The total transfer function T(s) for 
displacement becomes: - 


Tp (s)= cage (1.6) 
TIG a Zm) : 
m=1 

z, = (—.063,0.i) rad/s 22 = (—.063,0.i) rad/s 


24= (—4.443,4.443i) rad/s 24 = (-4.443,-4.443/) rad/s 


A = 1.23 *10' (Sum tz) 


Comparison with other broad-band systems, e.g. the Graefenberg array (Seidl & Stammler, 
1984), is simplified by showing the amplitude response to ground velocity in Vm7!/Hz as 


well as the phase response in degrees and group delay in seconds. The transfer function 
Ty(s) for velocity then becomes: 


Tp(s) 
Ds (1.7) 


Ty (s)= 


where D translates Volts to Counts (D= 104857.6 Counts/Volt). Figure 1.2.3 shows results. 


The main differences between the NARS and Graefenberg responses are compiled in table 
1.3 


Table 1.3 Comparison of broad-band systems 
Parameter Graefenberg NARS 
Gain 2400 Vs/m 3450 Vs/m 
High frequency -3dB point 5.0 Hz 1.0 Hz 
slope 42 dB/octave 12 dB/octave 
Low frequency -3dB point .05 Hz 01 Hz 
slope 12 dB/octave 12 dB/octave 


NARS is in general more sensitive to the lower frequencies and has only a 2™ order anti- 
alias filter compared to a 7” order for the Graefenberg instruments. Seidl & Stammler 
(1984) show that the anti-aliasing filter determines mainly the leading edge of the onset of a 
phase on the seismogram. The use of a 7” order anti-aliasing filter may introduce an 
apparent delay in the arrival time of 0.1 second at a noise level of 10%. This can be a seri- 
ous drawback, especially for a mobile array like NARS. In practice a 2” order anti-aliasing 
filter has proved to be sufficient for NARS, and results in an apparent delay well below the 
sampling interval even for high noise levels. 
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Figure 1.2.3 Individual instrument responses and total system response, 
showing: a. Amplitude of the transfer functions in Vs/m as a function of 
frequency. b. Phase shift in degrees as a function of frequency. c. Group 
delay in seconds as a function of frequency. The legenda shown in a. is also 


valid for b. 


1.3 SYSTEM CALIBRATION 


Calibration of the NARS stations is effected at every change of tape (roughly every 
14 days). A blockpulse of -1.0 Volt, with a duration of 32 seconds, is applied to the calibra- 
tion coil of the seismometers. Their response to this pulse, together with the pulse itself are 
recorded on cassette tape. An example of a recorded system calibration is given in figure 


(1.3.1) 


NEO4 Z comp 
FEB 43 (0443, 1986 
OL OF7-BF-E: 


NEQ4 NS comp 
FEB i3 (044), 


[Counts ] 


NEQ4 EW comp 
Fy 
09: 07: 27.216 


re) 
X 10+2 : 
time [s] 


Figure 1.3.1 Recording of a system calibration. Ticks in the EW component 
are minutemarks to check the timing. 


At system calibration the calibration coil applies a step in acceleration (Aa) to the 
system, followed after 32 seconds by a step (-Aa). This step in acceleration is given by: 


Aa = eat Feat (1.8) 
rv : 


Tai, the current through the calibration coil has a value of 4.5 mA. G,,; denotes the motor 
constant of the calibration coil with a value of 0.029+.002N/A and M stands for the weight 
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of the inertial mass of the seismometer, i.e. 2.0 kg. The system response to Aq is: 


Aa T; 
A(s)= 20 (1.9) 


and is given in Counts /Hz. 

There are a number of methods (see Rasson & de Meyer, 1983 for references) to 
measure the instrument parameters from the calibration pulse, so a check on instrument per- 
formance can be made. For the in this thesis presented higher-mode surface wave phase 
velocity measurements differences in magnification between subsequent stations are not 
important, since the spectra of each station used in the stacking procedure, see chapter 2, are 
normalized. However, for body wave modelling studies this is of crucial importance. We 
will show how one can calculate the magnification parameter A from the measured peak- 
to-peak amplitude of the calibration curve provided that one knows @p,@;,/) and h, from, 
e.g. , a Steady state calibration. Equation (1.9) can be written as: 


~—____S@A 
(S+@p)*[(s +@,h1)*+O7h7) 


A(s) Mp = .063Hz 50, = 6.283Hz 3h; =0.707 (1.10) 


In the time domain this gives (Jarosh and Curtis, 1973): 


2. 2 
a(t)= ea [2E +t (E74?)] e+ [-2Ecosnt + (S = ) sinne] em 
E=(@o-Mhy) N= hy (1.11) 


Since @)<«0,h, the term e~ is dominant after the first three samples for the NARS instru- 
mentation, thus we can make the approximation: 


AAa 


a(t)= Cand? (2E+t(E74?)] e (t>0.375s) (1.12) 
having a maximum at 

Gea Nies oe 

m= Oy (624?) 


For the NARS instrumentation the foregoing analysis has to be modified since the total cali- 
bration curve consists of two parts: A negative peak from the application of -1 Volt to the 
system, followed after t 32) seconds by a positive peak from switching off the calibration 
voltage. This sequence can be described by: 


a(t)=a(t)H(t) - a(t-t)H(t-1) (1.13) 


where H(t) denotes the Heaviside function. Differentiation with respect to t gives extrema 
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for this function. 


For 0< t <t: 
tin) = eae assuming tmini<t (see fig 1,3.2.); 
@o E741 
and for t2>t: 
bmax2 = tminit ie =Imini t+ OT 
1-e™* 


The peak to peak amplitude A,, is given by: 
App =| a (tmin1) ~ 4 (tmax2) | 


Gane 7 Doman [2EX + tnin (E24n)r +b 1(674n7) [A-1] 


—1(6747) e eb) | 


N=1-[l-e™ Je" (1.14) 


t is known accurately, t = 32.0s, and from ¢,,;,; and t,,..2 we can even check the parameter 
@ . Unfortunately the parameters ¢,,;,; and t,a,2 are hard to measure accurately since the 
calibration curve is very flat near its extrema. Microseismic noise disturbes the measure- 
ments considerably, especially in the stations in the Netherlands and Denmark. In practice 
the measured peak-to-peak amplitude of the calibration curve shows variations of up to 25% 
at the same station. This is due to station conditions, especially temperature and athmos- 
pheric pressure fluctuations. Figure 1.3.2 shows a calibration curve recorded in station 
NE06 (Witteveen) together with the theoretical response calculated from equation (1.11). 
Apart from a difference in absolute magnification the curves match well. Synthetic seismo- 
grams of fundamental mode Rayleigh waves, calculated for NARS instruments, show that 
the amplitudes are consistent with the recorded data within 20% (Snieder, pers. comm.), as 
was expected from the calibration results. A steady state calibration of all the stations is 
planned for early 1987. 
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Figure 1.3.2 Comparison between theoretical (T) and measured (M) 
calibration curve. (A) is the neglected factor [..] e~™ in equation (1.11). 


1.4 NARS INSTALLATION AND OPERATION 


As described in the general introduction, the NARS experiment has as its primary aim the 
measurement of the dispersion of higher-mode surface waves of the west-European plat- 
form. The higher-modes appear as an interference pattern in recorded seismograms and the 
dispersion of the individual modes can only be directly measured in the wavenumber 
domain through a stacking procedure (see chapter 2). For an optimum result it is necessary 
that all NARS stations are functioning at the occurence of an event. Therefore the installa- 
tion and operation of NARS described in this chapter is focussed on keeping as many sta- 
tions as possible functioning at the same time. 
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1.4.1 Organisation 


Thanks to the generous cooperation of institutes, station managers and country 
supervisors listed in table 1.4, we were able to carry out the NARS experiment. The task of 
a Station manager is to change a cassette tape once every 14 days and sent it to the data cen- 
tre in Utrecht. In case a station malfunctioning is detected the datacentre should be con- 
tacted as soon as possible, directly or with the aid of the country supervisor. 


Table 1.4 NARS station managers and country supervisors 

Station Station Manager 

G. Lind; 

V. A. Glargaard; 

J. Sorensen; 

L. Beuving; T. Jipping; P.J.M. Wit; 

A. van Wettum; 

J. Rasson; 

A. Simonin; 

J.-H. Pelletier; 

R. and B. Darchen; Mr. Bages 

M. Vadell; J. Esprabens; J. Bouillon; 

G. Cremades; M. Diez; 

M.A. and M.-J. Bordeje Muguerza; 

M. Herraiz; B. Benito; 

D. Hergueta; 

F, Vidal-Sanchez; F. de Miguel; G. Alguacil; J.M. Guirao and G Olivares 

F.L.M.V. Bergsteijn; 

J. Dorel; J. Fourvel; 

G. Payo; 

S. Dujancourt; 


Country Supervisor 


Denmark S. Gregersen; E. Jensen; 
France A. Simonin; N. Jobert 
Spain M. Herraiz; A. Udias 


1.4.2. Installation 


NARS began its operation in March 1982 with the installation of an experimental 
Station in Utrecht (NEO5). Within 4 months also the station in Dourbes (Belgium, NE06), 
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and Witteveen (NE04) were installed. With these three stations within a radius of 200 km 
from the datacentre in Utrecht we were able to gain some experience with the instrumenta- 
tion (and with possible customs problems in west European countries). 


Figure 1.4.1. Installation scheme of the NARS array. [A] June 1982, [B] 
November 1982, [C] January 1983, [D] May 1983, [E] November 1983 array 
complete. Open triangles denote stations that were closed during 1983-1986 
and moved to stations NE15-NE18, [F] December 1986. 


In November 1982 three stations could be placed in France, NE08-NE10, in January 1983 
two stations in Denmark and one in Sweden, NEO1-NE03, and in May 1983 four stations in 
Spain, NE11-NE14. Finally the array became complete with the installation of NEO7 in 
November 1983. Figure 1.4.1 shows the different stages of installation and table 1.5 the sta- 
tion locations and their period of operation. In six cases:NEO1,NE05,NE08,NE09,NEi1 
and NE12, station locations proved to be inadequate after some time and alternative loca- 
tions had to be found. 


17 


Table 1.5: NARS station coordinates. 


Coordinates Elevation(m] Period of operation 
NEO1 Gothenborg 57.801N 12.132E 55.0 Feb. 83- Jan. 86 
NE02 Monsted 56.459N 9.170E 60.0 Feb. 83 
NEO3 Logumkloster 55.045N 9.153E 25.0 Feb. 83 
NE04 Witteveen 52.813N 6.668E 17.0 Jul. 82 
NEOS Utrecht 52.088N 5.172E 2.0 Mar. 82! 
NE06 Dourbes 50.097N 4.595E 225.0 Jul. 82 
NEO7 Villiers-Adam 49.074N 2.232E 70.0 Nov. 83 
NEO8 Aigurande 46.420N 1.730E 360.0 Nov. 82- Dec. 84 
NE09 Les Eyzies 44.852N 0.981E 160.0 Nov. 82- Feb. 86 
NE10 Arette 43,086N 0.699W 480.0 Nov. 82 
NE11 La Almunia 41.477N 1.372W 370.0 May 83- Nov. 83 
NE11 Ainzon 41.814N 1.517W 440.0 Nov. 83 
NE12 Valle delosCaidos 40.642N 4.155W 1280.0 May 83- Feb. 85 
NE13 Puertollano 38.685N 4.091W 700.0 May 83 
NE14 Granada 37.190N 3.595W 774.0 May 83 
NE15 Valkenburg 50.867N 5.785E 100.0 Jun. 84 
NE16 Clermont Ferrand 45.763N 3.103E 80.0 Nov. 84 
NE17 Toledo 39.881N 4.049W 480.0 Feb. 85 
NE18 Rejaudoux 45.304N 1.516E 410.0 Feb. 86 


1 Station closed between June 1984 and Januari 1986 


1.4.3. Data processing 


Until May 1985 NARS data have been processed on a HP-1000 mini computer with 


limited disk space (120 Mbyte) and one 1600 bpi tape unit. The NARS data format had 
been chosen to satisfy our needs within these constraints, but was incompatible with other 
existing data formats. 

In May 1985 the HP-1000 was exchanged for a GOULD PN6080, a 32 bit machine 
equipped with a UNIX Berkely 4.2 operating system, 1020 Mbyte disk space and two 
1600/6250 bpi tape units. 

The data format in use since then is compatible with the Global Digital Seismograph Net- 
work (GDSN) format (Hoffman, 1980). This facilitates the exchange of data and enables us 
to use all data distributed by the GDSN together with NARS data, without any format 
conversions. Furthermore, a vast amount of software developed by the GDSN for retrieval 
of the data could be used. Because NARS is the first large scale broad-band experiment in 
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global seismology, we will describe the data processing in detail. 

After the cassette tapes with the data are received by mail, a rough print-out of the data is 
made with the aid of the play-back unit (CTU-1). Obvious instrumental problems are 
immeadiately detected and action can be taken for repair. Cassette tapes are subsequently 
read on a Micro-Sol microcomputer equipped with an CP/M operating system and C- 
compiler. Data are stored on a 20 Mbyte Winchester disk and some preliminary calcula- 
tions are made. These calculations provide insight in the nature of the recorded triggering. 
Routinely a check is made on: 


[1] The calibration, peak-to-peak amplitude is calculated. 
(2] Microseismic noise, dominant frequency is calculated 
[3] Parasitic signals, dominant frequency is calculated. 
[4] Identification of the seismic event. 


For each triggering start and end times, instrument status information and parameter settings 
are written to a report file together with outcome of the preliminary calculations mentioned 
above. Report files are printed and kept in a binder. Valid data are transmitted from the 
MicroSol to the Gould system and converted to GDSN format. This means that the 114 
byte PDR-2 header is converted to a 20 byte GDSN header. No information is lost, since 
most of the PDR-2 header contains status information on the PDR-2 system and this infor- 
mation is kept in the report files. Data remain in 16 bit multiplexed format. The mantisse is 
converted to a two’s complement code and the exponent changed from a maximum of 2’ to 
2!°. After conversion data are stored on disk in station directories and plots of the first 400 
seconds of all accepted events are made and kept in binders. When data from all stations of 
a certain time period are processed event tapes are constructed. Apart from the data, pre- 
ceded by a data-log that contains information on system response and calibration, a station- 
log and event-log are constructed. 

The station-log contains not only the location of the station, but also the type of data (inter- 
mediate period) and time corrections to be applied to the data. Time corrections are rou- 
tinely determined by measuring the onsets of the minute pulses on the EW channel. 

The event-log gives a list of station names that recorded this event together with the position 
of the recordings at the tape. Selection criteria for events to be stored on event tapes are: 


{1] Events with magnitude 25.8 
[2] Events of magnitude <5.8, but recorded at 3 or more stations. 


[3] All European events that are reported by the Preliminary Determination of Epicenters 
(PDE), compiled by the U.S. Geological Survey. 


These criteria are in use from 1985 onward. Event tapes from 1983 to 1985 are assembled 
without the third criterium and with the treshold magnitude for the first and second criterium 
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set to 6.0. Apart from using the event tapes in our own laboratory, event tapes are shipped 
to the USGS, where the data are incorporated in the GDSN data base and distributed. 

For NARS data a separate disk of 300 Mbytes is available at the Gould system. In practice 
approximately 80 Mbyte is used for PDR-2 data handling and 200 Mbytes can be used for 
data analysis. Selected events can be retrieved from tape and are permanently available to 
users. For graphical display and simple data manipulation of the data the Seismic Analysis 
Code (SAC) developed by the Lawrence Livermore Laboratories is in use. Data retrieved 
from event tape are converted to the SAC format for this purpose. 

For handling of file errors (such as polarity reversals or misorientation of horizontal com- 
ponents) a patch file is maintained that contains all information needed to correct the data 
when being converted to SAC format. This way, the original data (on tape) are always left 
untouched. 


1.4.4. NARS Performance 


The performance of the NARS array is measured by selection of events with a body wave 
magnitude higher than a specified treshold and counting the number of stations that 
recorded this event. The use of event logs greatly facilitates this search procedure. Figure 
1.4.2 shows results of this procedure. 

In the installing phase our policy was to carry out station repair as soon as possible. This 
resulted in a high average return rate with a peak at approximately 400 days (early 1984). 
The return rate is defined as the number of stations that succesfully recorded one specific 
event. In this period station repair was too time consuming and we changed to a system 
where station servicing was limited to two or three times a year. In this way all stations 
were expected to work optimally at least one to two months after servicing and most contin- 
ued doing so for a larger time span. This change resulted in a lower return rate, but still a 
considerable number of events were recorded in 7 or more stations. This is a minimum 
requirement for higher mode surface wave analysis using NARS. 

In the new set up, as explained in the last chapter, European events are incorporated even if 
they are recorded in only one station. This is the reason why a concentration of mainly low 
magnitude events recorded in less than 3 stations is visible after 800 days. A catalogue of 
all recorded events presently available on event tape is presented in appendix C. 
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Figure 14.2. NARS performance. 
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blackness indicate events with body wave magnitude [1] >4.0; [2] >6.0; [3] 
>6.4, The black line indicates the number of stations that were installed at the 


time of measurement. 


Station repair has been necessary because of many instrument failures. A description of the 
major problems is given below. 


Event recorder, 


Software problems: 


During the initial stage of the project (3 stations running) two problems were detected: 


(1) 


The original idea that the triggering algorithm should work on all three components 


simultaneously turned out to produce too many false triggers. It was decided to 
change the trigger to work on the vertical seismometer only. 
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[2] STA/LTA ratio triggering is sensitive to a DC offset of the pre-amplifier output. A 
modification of the algorithm, consisting of subtraction of the true LTA before taking 
the absolute value of the samples to calculate the absolute LTA, described by Proth- 
ero (1980), can solve the problem in a next generation event recorders. 


Hardware problems: 


During 1982,1983 and 1984 station operation was hindered by many failures in micro- 
electronic components. Especially EPROMs (Erasable Programmable Read Only Memory) 
and data RAMs (Random Access Memory) broke down regularly. In the first 500 days of 
operation at least 8 EPROMs and 10 data RAMs had to be replaced. An additional source 
of trouble is the use of sockets to mount the EPROMs and data RAMs. Pins incorrectly 
mounted occasionally make contact and generate errors. Sockets were removed for the data 
RAMs in 12 PDR-2’s. During 1985 and 1986 the EPROM problems diminished, but RAM 
defects kept troubling us. 

Installation of stations in a very humid surrounding, like a cave or old mine, causes corro- 
sion troubles. Some parts of the tape-drive of the event recorder were originally not plated. 
Replacement with plated parts in subsequently all stations solved this problem. 

Information is recorded on cassette tape using four tracks, three signal and one parity track. 
Early 1986 wires connected to the write-head of the recorder came loose in two stations. 
When reading these tapes at the data centre the play- back unit showed good quality data for 
two of the three channels, the third one showing no data at all. Data transfer to the micro 
computer however was not possible, the recorder could not find any. data. This problem 
could be solved by sharing the fourth track with the absent channel, so the recorder could 
synchronize at the start of a cassette tape. Since the header is (partly) unreadable in this pro- 
cedure, the recording date and time is not present. If the absent channel is the EW com- 
ponent, time can only be roughly reconstructed by lining up of onsets from different sta- 
tions. 


Seismometers, 


The sudden appearance of parasites in the seismic signal, in several cases after two years of 
operation, did upset our triggering and resulted in filling one cassette in one day. These 
parasites only occur on the vertical channel. An example of such a feature can be seen in 
figure 1.4.3. The fact that only the vertical instrument was affected limited the problem to 
the spring or flexures. This problem is under study at the manufacturer. 

Humidity caused oxidation problems mainly for the horizontal seismometers. These instru- 
ments are not sealed with a rubber ring as the vertical component. The casing of the instru- 
ments was modified as that of the vertical seismometer and the horizontal instruments were 
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replaced. An additional effect of humidity on seismometers is the generation of an offset. 
The connection between the seismometer suspension and the casing consists of very thin 
copper wires mounted on a terminal strip. Condensation of water particles on this strip 
results in leakage between the three components. Humidity problems could also be dimin- 
ished by surrounding the seismometers and pre-amplifier with a blanket made of very thin 
aluminum foil. The small amount of heat produced by the power supply of the pre-amplifier 
prevents condensation of water particles on the seismometers. 
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Figure 14.3. An example of a long period parasitic signal, only visible on 
the vertical component. 


The sensitivity of the vertical seismometer to temperature changes is a major hindrance in 
station operation. A change of only 3°C is sufficient to move the mass so far out of its 
equilibrium that proper functioning is impossible. The change in outside temperature from 
summer to winter and vice versa can easily be 30°C along the array. Since the array should 
be portable and therefore not restricted to well isolated environments we had to find a solu- 
tion that can be installed and removed in one day. 

For this purpose a temperature isolation box has been developed that consists of two polys- 
tyrene covers with a small heating cell inside the top of the outer cover in order to minimize 
turbulent air currents (see Willmore, 1979). Direct radiation from the heating cell some- 
times disturbed the seismic signal. This problem has been overcome by using tiles glued to 
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the outside of the inner box. A detailed description of the box is given by van der Meer (in 
prep.). 


1.5 CONCLUSION 


From our experience in deploying a mobile, digital, broad-band seismograph net- 
work like NARS we conclude that it is feasible to carry out a seismic experiment with only 
a very limited staff. The help of station managers and country supervisors in maintaining 
the stations and establishing a good communication between the datacentre and the seismic 
stations is of vital importance. 

One drawback of the very rapid development of (digital) instrumentation is the fact 

that the products seem to be not very well tested. The continuous problems we encountered 
in failing components like EPROMs and data RAMs and some construction errors support 
this conclusion. As a consequence we think it very important to be able to monitor the 
status of a remote station directly from the data centre. This will be a future development 
for the NARS instrumentation. 
For this study the consequences of the instrument failures are that the original idea of a 
sufficient station coverage is not fulfilled. Therefore it is advisable to look for algorithms 
that can enhance resolution in the process of stacking of higher mode surface waves using 
sparse datasets. This will be explored in the next chapter. 
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PART 2: ANALYSIS OF HIGHER MODES OF SURFACE WAVES 


As was discussed in the general introduction, the measurement of higher-mode sur- 
face wave dispersion became possible with the introduction of a stacking technique (Nolet, 
1975). A prerequisite for the application of this technique is the existence of a linear array 
of seismographs like NARS. In an analysis of the possibilities and limitations of this tech- 
nique, Nolet and Panza (1976) came to the conclusion that for regions of mild lateral hetero- 
geneity the main source of error in the measurement of phase velocities was due to the 
inadequacy of the array response. Since the array response is fully determined by the posi- 
tion of the stations, NARS has been designed as to minimize the influence of an inadequate 
array response. In the design a return rate of 80% was expected for NARS. However, in 
practice the average return rate of the data did not exceed 60-70% (see chapter 1.5.4). 
Therefore we decided to investigate methods to enhance resolution. 


2.1 RESOLUTION ENHANCEMENT USING A PRIORI KNOWLEDGE 


Higher mode surface waves show up in the seismogram as a pattern of interfering 
waves arriving in the time interval between the S onset and the fundamental mode surface 
wave, In order to separate the different modes we have to transform the frequency spec- 
trum w (A,q@) of the signal s(A,t) to the complex wavenumber spectrum W (k ,@): 


Cy 


W (ko) = | w(Aw)e*4 dA (2.1) 
where & denotes the wavenumber, @ angular frequency and A epicentral distance. By 
showing the equivalence between surface waves and normal modes in the limit for large 
radial orders, Nolet (1976) derived an expression for the frequency spectrum: 


w(Aj.0) =F Fyj(a) ero (2.2) 
n=l 

F,,;(@) is the amplitude of mode n in station j and k, (@) the wavenumber of mode n. @,; is 
the initial phase and depends on the azimuth from the source. The most important approxi- 
mation made in the derivation of equation (2.2) is the asymptotic representation of a Legen- 
dre polynominal for A>>0. This approximation imposes limits to the location of the array 
with respect to the source. Panza et al (1973) determined as condition for a 3 figure accu- 
racy in amplitude and phases: k A>10, which is appropriate for the measurements presented 
in this thesis. If stations do not differ considerably in azimuth and if the array is not situated 
near a minimum of the radiation pattern we can neglect the variations of 6,,;(@) in the dif- 
ferent stations (Nolet and Panza, 1976). Amplitude F,,; is a function of station position, but 
the effect of dissipation over the array is small. We will therefore neglect the station 
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dependence in the amplitude. If necessary, we can correct for the geometrical spreading. 
Combination of equation (2.1) and equation (2.2) without station dependence of amplitude 
and initial phase gives: 


M : 
W (k,@) = 20 ¥ F,(@) e'* © 8k, Kk) (2.3) 
n=l 
The complex wavenumber spectrum consists of a set of M delta functions 5(k, —k ), centered 
around the modes n. Unfortunately, , we have a finite amount of recording stations and there- 
fore we can only make an estimate W (k) of W(k): 


2 1X F 
Wk.o) = = za (Aj,) w(Aj,0) eo (2.4) 
jz 
The contribution of each station j is weighted by a factor a (A; ,@). Combination of equation 
(2.4) and (2.2) gives at constant frequency: 


a M ‘ 
W(k)= > F, e'* H(k,-k) 


n=l 


H(k)= pe y pikAy 

N p2 a(Aj)e (2.5) 

j=l 

|H(k)\? is called the array response function and its shape determines the resolution in the 
measurements. Because the number of stations is very limited, H(k) does not resemble a 
delta function as in (2.3)- which we try to approximate. Although |#(k)I is peaked in 
k=0, there are many secondary maxima. This is essentially an effect of spatial aliasing and 
we Shall refer to these peaks as ‘aliasing lobes’. See figure 2.1.1 for a typical example of an 
array response function shape. Since we know the location of the stations, and thus the 
shape of the array response function, we can correct the wavenumber spectrum for the 
effect of spatial aliasing and possible non-equidistant sampling. The CLEAN algorithm, 
introduced in seismology by Nolet and Panza (1976), iteratively carries out this procedure. 
Cara (1978) shows the usefulness of choosing weights [a(A;)] in order to manipulate the 
array response function. To achieve this he designed a penalty function that was subse- 
quently minimized. Accurate determination of the phase velocity of mode n can be 
obtained if the number of modes present in the signal is small. Therefore small segments of 
the seismic signal are selected in the time domain around a traveltime t = A,/v, and estima- 
tor (2.4) calculated at a chosen frequency (Wp ). This procedure is repeated for a suite of 
group velocities and results in a mode separation diagram of phase- versus group- velocity 
for Wo (see figure 2.3.8). 
Our present aim is to investigate whether the above described stacking technique can be 
improved by using optimum filtering techniques to determine the weights a(A,) or by 
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modifying the estimator equation (2.4). 
2.1.1. Optimum weights 


A comparison with classical deconvolution provides insight in the role of the weights a (A;) 
in equation (2.5). For this purpose we will introduce a stationary time series a(t) that can 
be written as a convolution of a wavelet B (t) with a series of spikes C, : 


M 
a(t)= 5 C, B(t-t,) 


n=} 


M 

UGS, F, G(k-k,), G(k)=H"(k), F,=F, e'* (2.6) 

n= 

From equation (2.6) one can observe a strong analogy of the estimator W(k) with a(t). As 
a consequence our problem of recovering the modes from our estimator is equivalent with 
wavelet deconvolution. The complex conjugate of the array response function [G (kK )] is 
analogous to the wavelet. This allows us to apply a vast amount of well known filtering 
techniques (e.g. Robinson and Treitel, 1980; Oldenburg, 1981). One important difference 
we should keep in mind in this comparison is that the station location of NARS is not 
equidistant. In Appendix B we will show the importance of choosing non-equidistant pat- 
terns and how an optimum array configuration can be determined. 
The role of G (kK) can be clarified with a well known excercise from Fourier theory. If we 
define the discrete Fourier transform of the finite series x (A;),j=1,.....N as: 


X(k)= a x(Aje" (2.7) 


x(A;)= J X(k) elt 


(2.8) 
it can be shown by substitution of (2.8) in (2.7), that the true X(k) can not be uniquely 
reconstructed from the estimated X (kK), due to the (un)weighted array response function or 
its complex conjugate G (kK): 


Xk)= [XEIGE-*) a (2.9) 


We assume that X(k) is bandlimited (-K<k <x). This assumption is not unreasonable in 
seismological practice and serves to avoid or reduce problems of spatial aliasing. The shape 
of the array response function (or its complex conjugate) is essential in our reconstruction 
and we can try to manipulate it by minimizing a form like: 
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K 
JIG &)}-LE)I ak (2.10) 
—-* 
where L(k) is the target function we want G(k) to resemble. Usually, this function is 
chosen as a delta function or a Gaussian function with prescribed width. 
Before we treat this problem in detail it is important to understand the relation between the 
choice of weights in the A domain and a filter operation in the wavenumber domain. It will 
be shown explicitely that these two operations are equivalent. In the rest of this chapter we 
will write functions in the wavenumber domain with a capital, their representation in the A- 
domain in lower case. 
If we denote the unweighted ’wavelet’ with G’ (k): 


; 1 oe 
G k= he : (2.11) 
j-l 


we can show that the application of weights a(A;) is a filtering operation in the k-domain 
and can be written as a convolution with a filter response V (k): 


G(k)=G'(k) * Vik) (2.12) 


* J G’ (K) V(k-) dK 


Substitution of (2.11) and using (2.7) leads to: 


spe ~ikAn 
Gtk)= W Dd van) e 2n 5(A,,—A;) 


m=) 


N 
= a Lv aje"™ (2.13) 
j=l 

And we see that the factor 2xv(A;) can be identified with the weights a(A;) in equation 
(2.5). We will now follow the theory of Byrne and Fitzgerald (1982,1984) to show the 
influence of the band-limitness of a signal in an elegant way. A function P (k) is intro- 
duced, which contains a priori known information of the spectral density of the signal. 
Furthermore we consider the Hilbert space H = L?(—ce,00;1/P (k)) with inner product 


< ¥(k),Z(k) >» = f YO)ZH aa x (2.14) 


~oo 


Suppose G (k) is band-limited on [—«;«], then we can minimize the weighted norm: 
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t 21 a z 
J P®)6®)-L@)! PH dk P(k)=1+e for -KSk <K 


=e for k elsewhere (2.15) 


€ is introduced to incorporate an average level outside the limited bandwidth. This parame- 
ter plays the role of a damping parameter. After minimization of equation (2.15) with 
respect to a” (A,,) we obtain: 


N 
D 4@(An) P(A, — Am) =N 1(A,) (2.16) 
m=1 
with p (A;) and /(A;) the inverse Fourier transform of P (k) and L(k) respectively. For our 
choice of P (k) we can write: 


sink (A, — A,,) 


pA, ee) aoa ee —A,,) (2.17) 


In case of equidistant sampling, A; =j A, and a band-limitness that coincides with the 
Nyquist interval, k= 2 1 / A, equation (2.16) simplifies to: 


a= oN 1(A,) (2.18) 


The highest resolution is obtained when we take for L(k) a delta function, so /(A,) = 1 and 
weights are constant. This means that for equidistant location we can not expect an increase 
in resolution by using weights. In any other case we have to solve equation (2.16) to calcu- 
late weights a(A;). For L(k) = 8(k) we will further investigate the influence of the weights 
a(A;) on the shape of the array response. Of main interest are two parameters: 


{1] The width of the main lobe of the array response (determines the resolution). 


[2] The height of the aliasing lobes. Aliasing lobes dominate the ‘noise level’ in the k- 
domain 


As can be seen from equation (2.17), the parameters that determine the weights are the 
band-width « and damping parameter ¢. In figure 2.1.1 we show results for the 14 station 
NARS configuration. In this figure the width of the mainlobe and height of the first side 
lobe is measured as a function of « for 3 values of €. Two regions can be observed: 

[a] «>1.2107: 


In this range the damping parameter has no influence on the shape of the array response 
function. A trade-off exists between the width of the main lobe and the height of the first 
side lobe. 
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[b] «<12107: 


The width of the mainlobe decreases with decreasing «. Increase of the damping factor 
slows down this phenomenon. The height of the first sidelobe oscillates for ¢ = 0 as a func- 
tion of « and behaves much smoother when damping is increased. When « becomes smaller 
than the width of the main lobe (approximately given by k,), the algorithm becomes 
unstable. 

The region k,<«<1.2 10~ is favourable for increasing the resolution without blowing up the 
first side lobe. The question rises what happens to the other side lobes. To investigate this 
we choose a damping factor of 1% and calculate the height of the first 5 side lobes as a 
function of «x. Results are shown in figure 2.1.2. 
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Figure 2.1.1. The top figure gives the unweighted array response function |G (k) 17. 
Indicated are the Nyquist wavenumber k, = 27/A and the width of the mainlobe 
Kk. = 2n/NA, where A is the average station distance. Marked with a 0 is the height of 
the first minimum aside the main lobe, and 1 the height of the first side lobe. In the 
lower figure the solid lines give the position of the first minimum (k-width) as a function 
of band width « for: [1] €=0% ; [2] €=.2% and [3] €= 1%. The percentages are 
relative to p (0). Dashed lines give the height of the first side lobe (G-max) as a function 
of «, at the same values of € as the solid lines. Horizontal lines indicate the values for 
the unweighted situation and act as reference lines. Labels 0 and 1 in the upper figure 


correspond with the labels at the horizontal lines in the lower figure 
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Figure 2.1.2. The description of this figure is the same as figure 2.1.1., except that the 
width of the mainlobe is omitted and that the dashed lines now give the height of the first 
5 lobes. The damping parameter € = 1%. 


The region 1.2 10? < «x $2.4 10° is favourable to diminish the width of the first five lobes. 
Combination of the outcome of these numerical experiments leads to the following conclu- 
sions: 
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[1] The width of the mainlobe can be significantly diminished for small values of «. 
However, the height of the subsequent sidelobes increases rapidly with decreasing x. 


[2] A region for K can be given where the first five sidelobes are significantly diminished 
in height with respect to the unweighted array response, in the case of figure 2.1.2. 
this can be realized by taking 1.2 10? <« <2.4 107. 


Constructions of figures like 1.2.1 and 1,2.2 for a specific station distribution can be helpful 
to design an array response function. The station distribution will be the limiting factor in 
the usefulness of applying weights. 


2.1.2. Fourier estimators incorporating a priori knowledge 


As we have seen in the previous chapter, the incorporation of a priori knowledge 
provides us with a powerful tool to enhance resolution. In this chapter we will show how we 
can use a priori knowledge about the spectrum in the Fourier estimator itself. This theory 
has recently been developed by Byme and Fitzgerald (1982,1984) and Michette et al (1984). 
Since we know the approximate location of our modes from previous studies, these tech- 
niques are expected to improve the accuracy of the measurements. Suppose we put our a 
priori knowledge about the shape of the wavenumber spectrum W (x) in a function P (k): 


P(k)=A,+€ kyz-K, Sk Sk,+K, , n=0,1,2.... 


=E elsewhere (2.19) 


with k, the a priori known location of mode 7, «,, the uncertainty in the a priori knowledge, 
A,, the amplitude of mode n and e a damping parameter (see last chapter and figure 2.1.3). 
If we assume W (k) to be a member of the weighted Hilbertspace H = L4Rs090; 57) with 
inner product (2.14) we can write: 


= <W(k).8;(k)>> 


, ina, ak 
w(aj)= J Wkde RE eon 


™ 


gj(k)=e *” P(k) (2.20) 
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Figure 2.1.3. The absolute value of the complex wavenumber spectrum for one row of a 
mode separation diagram. The vertical scale is in dB. The maximum of the entire 
diagram has been set to 100. The dashed curve denotes the conventional DFT stack, the 
solid curve the PDFT stack. Superposed is the function P(X) in arbitrary units. No 
information about relative excitation of the different modes is incorporated, i.e. boxes 
have equal height. Arrows mark the exact locations of the modes. 


An estimator W (k) can now be constructed, as in Backus Gilbert theory, using the data ker- 
nels g, as a set of basisfunctions: 


— N 
Wk) = 1 558;) (2.21) 
j=l 


In contrast to the estimator (2.4), based on a Discrete Fourier Transform (DFT), our param- 
eters b; are not equivalent to the (un)weighted spectral values w(A,). If we want to con- 
struct our estimator from the data w (A; ) we have to minimize the weighted norm: 
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=e 


Pb) dk (2.22) 


fiw) - DojpP&e 8" 1? 
eo j=l 


Differentiation with respect to b, leads us to 


* | re (ana) ak 
Wkje*™™ = yb, [Pek awY & 
J wane Pe hae ie 
or 
N 
j=l 


Coefficients b; can be determined. Our choice for P(k) leads to equation (2.24) for the 
function p (A, —A;) in (2.23) 


M it. ta,-a.) Sin K(A;—A,) . 
p(,-A;) = >» A, eit a) aay e6(A; — A,) (2.24) 
m= j n 
and the estimator becomes: 
= is -ik A, 
W(k)=P(k) Hoje” (2.25) 


j=l 

This estimator is called the PDFT estimator (Byme and Fitzgerald, 1984) and is a generali- 
zation of the band-limited extrapolation procedures (Papoulis, 1975). 

The stacking operation now consists of two steps: first one has to calculate the coefficients 
b, from equations (2.23) and (2.24), then these coefficients are summed with their appropri- 
ate kernels g;(k) in equation (2.25). Figure 2.1.3 shows an example of the enhancement of 
energy in one row of a mode separation diagram when we use equation (2.25) in stead of 
equation (2.4). 
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2.1.3. Comparison of algorithm performance 


A comparison of the different algorithms developed in the previous chapters with the 
conventional stacking technique can be made by measuring their performance on a set of 
higher mode Rayleigh wave synthetics. Since we know the exact values of the phase velo- 
cities for the different modes as a function of frequency, we can construct mode separation 
diagrams, measure phase velocities for the modes involved, and determine the error in the 
measurements. For these calculations we calculated synthetic seismograms using parame- 
ters identical to the March 24, 1984 Kurile event (see table 3.1) that was actually recorded 
in 10 NARS stations. The M7 model (Nolet, 1977) and centroid moment tensor solution 
published in the PDE monthly bulletins were used in the computation. 


2.52 1.56 -1.43 
M=| 1.56 4.82 1.96] 10'9 [Nm] 
-1.43 1.96 3.39 


The synthetics consist of fundamental and first seven higher modes, calculated in the fre- 
quency band between 10 and 100 mHz and are shown in figure 2.3.1. Before we can pro- 
duce mode separation diagrams for this comparison we have to investigate the possibility of 
using weights and the influence of a priori knowledge on the PDFT estimator. 


Weights: 


The usefulness of weighting can readily be explored by construction of weight diagrams 
like 2.1.1 and 2.1.2 for the station configuration of the test signal. Figure 2.3.2 shows that 
the first aliasing lobe can be significantly reduced in size for ky < Km < 2.4 10°, while the 
mainlobe width can only be reduced for «,, close to k2. The other sidelobes, however (see 
figure 2.3.3) blow up very quickly and there is no region for K where we can significantly 
improve the average height of the first five lobes together. These results are based on a 
configuration for only one event, but we selected the event/station configuration with the 
best coverage and the lowest average sidelobe level among all other configurations studied 
in this thesis. Construction of weight diagrams for the other events in table 3.1 shows that 
the discussed results represent the best case of all presently available events recorded by the 
NARS array. In our frequency range of interest, roughly between 20 and 70 mHz, the 
difference in wavenumber and group velocity between succesive modes is such that the 
effect of blowing up sidelobes can easily influence the measurements. In general, group 
velocity curves of the various modes approach each other as the frequency increases up to 
approximately 0.1 Hz. Cara’s (1978) measurements should also suffer from this effect at 
higher frequencies. The mode separation diagram he shows as an example for his choice of 
weights is calculated at a frequency of 14 mHz (70 s). The second higher mode that 
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epicentral distance |km | 


becomes visible after weighting is clearly separated in group velocity from the much better 
excited first higher mode. For his example weighting is an excellent tool, but application of 
the method to higher frequencies is rather doubtful. We can conclude that assigning 
weights to the station does not provide a useful tool in our case where we want to calculate 
many mode separation diagrams on a routine basis. However, for special purposes, e.g. to 
look in the vicinity of a specific mode, weights can be optimally chosen. 
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Figure 2.3.1. Synthetic Rayleigh wave data set. Parameters-from March, 24 1984 


Kuriles event (PDE). 
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Figure 2.3.2. The description of this figure is the same as in figure 2.1.1. 
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Figure 2.3.3. The description of this figure is the same as in figure 2.1.2. 


Influence of a priori knowledge on the PDFT estimator: 


When using a priori information in reconstruction of a spectrum by means of the PDFT esti- 
mator we have to be sure that this does not introduce a bias in the reconstruction itself. In 
order to check this, a test has been designed. The set up of the test is illustrated by figure 
2.3.4. 
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Figure 2.3.4. Test set up. Open arrows with solid shafts are true modes, at locations 
k, ,» Solid arrows with solid shafts are a priori assumed modes, at k, ,, and open arrows 
with dashed shafts are measured modes at k,, ,. For sake of clarity only two measured 
modes are shown. Parameters of interest are dk», the difference between true and a 
priori location and dk, the difference between true and measured location. Uncertainty 
in the a priori knowledge (box width) is given by Ak, ,. Modes are counted from right 
to left as phase velocity increases that way. 


The a priori location k,,, of the modes is in this experiment assumed to be located at a dis- 
tance (-1)”Ak,,, from the true location of the modes. This gives a direct relationship 
between k,, and k, ,: 


Kain 7 kin + (-1)” Ak, (2.26) 


P (k) is constructed as a superposition of boxcar functions around the assumed locations, 
while no amplitude information has been involved, i.e. all boxes have equal height or A,,= 
constant in equation (2.19). Parameters that determine the shape of P(K) are the uncer- 
tainty in the a priori knowledge, Ak, ,, and the average level ¢. € is specified relative to the 
height of the boxes. With these parameters we can compute mode separation diagrams and 
plot the difference between true and assumed mode location: 


dk» = Kran ie k, , (-1)"*1Ak, (2.27) 


for each mode n versus the difference between true and measured mode location: 
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dk =kin —knn (2.28) 


If dk, =dk», then the measurement is entirely determined by the a priori assumption 
(kan =kmn)- On the other hand dk, = 0 means that the a priori knowledge did not influence 
the result at all. In our application we want Idk,| to be close to zero, the measurement 
should be close to the true value, and |dk2| as large as possible to accomodate errors in the 
a priori knowledge. The design of our test is in the wavenumber domain, but the informa- 
tion we are interested in is how well phase velocities can be measured. Phase velocity vy is 


WO ¢,. P : persed 
related to wavenumber by vy = Since w is constant we can write for a perturbation in v, : 


yy =e (2.29) 
Relative changes in phase velocity and wavenumber are equal (apart from a sign). For this 
reason the parameters Ak, , and Ak, , are specified in percentage relative to k,, and k, , 
respectively. 

Three examples of such diagrams for a perturbation Ak, ,= 1% , damping e= 1% and dif- 
ferent values of “box width’ Ak, ,, are shown in figure 2.3.5. dk, and dk, are transformed to 
phase velocity values in [km/s]. This facilitates the interpretation of the diagrams in terms 
of phase velocity differences. Two clusters of points are visible, one for modes with a posi- 
tive and one for modes with a negative offset. The mean and variance of these two clusters 
are calculated and plotted as a function of box width Ak, , in figure 2.3.6a and as a function 
of damping (in a calculation with Ak, ,= 1% and Ak, ,= 2%) in figure 2.3.6b. In the follow- 
ing we will discuss the influence of a priori knowledge on the measurements with the aid of 
figures 2.3.5 and 2.3.6. 

If Ak, ,<Ak,, (the box width of P (x) is less than the difference between true and assumed 
mode location) the a priori location determines the final solution, but this influence dimin- 
ishes rapidly with increasing box width. For a box width of 2% is dk, on average close to 
zero. For Ak, ,=2% the results remain essentially the same. In realistic situations the error 
in higher mode phase velocities is assumed to be 1 to 2% (Nolet and Panza, 1976). We 
therefore adopt a conservative box width of 3% or more. This choice minimizes the 
influence of a priori knowledge. 

The dependence of ¢ is illustrated by figure 2.3.6b. A minimum is detected near e= 0.1%, 
and from this value the error dk, increases rapidly with decreasing value of . The algo- 
rithm becomes unstable. If e> 0.1% performance is somewhat worse, but remains within 
reasonable limits (.02 km/s). 
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Figure 2.3.5. Test results for Ak, ,=1%, damping €=1% and box widths Ak, , =0.5%, 
1% and 2%. The lower three figures show dk» as a function of dk,, the upper three 
figures show the situation around one specific mode. The dashed arrow represents the 
measurement. 


One can conclude that a priori information can be incorporated in the stacking pro- 


cedure without biasing the result. The only constraint is that the uncertainty in the a priori 
knowledge, i.e. the width of the boxes in the a priori function, is larger than the expected 
error in the measurements and that the damping parameter € is not too small. 
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Figure 2.3.6. Difference between true and measured mode location (dk ) averaged over 
measurements of 7 modes at 14 frequencies as a function of box width dk, , (top figure) 
and as a function of damping parameter € (bottom figure). Vertical bars indicate the 
standard deviation. Values for dk 2 are restricted to two regions at +[.052,.054]. 
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The choice of the damping parameter depends on the amount of noise present in the signal. 
The main source of noise comes from the aliasing lobes of the array response function and 
is dependent on station configuration, so no unique value can be assigned to €, its optimum 
value should be experimentally determined in each analysis. In the next section we will 
compare different algorithms and the effect of the PDFT estimator on the mode separation 
diagram will be discussed. 


Comparison of proposed algorithms: 


Now we are in the situation where we can compare the performance of the different algo- 
rithms on the synthetic dataset. This comparison can be carried out in the same way as 
described in the last section, namely by measuring the difference between true and meas- 
ured phase velocities: Avy. The algorithms we will use are: 

[a] DFT stack 

[b] DFT stack with CLEAN 

[c] PDFT stack 


One mode separation diagram, calculated for a constant frequency value, produces m phase 
velocity values for m modes. The mean and variance of these values are plotted in figure 
2.3.7 as a function of frequency, or in this case: period. From this figure we can see that the 
PDFT and original DFT stack give comparable results, but that CLEAN produces outliers 
for periods lower than 30 seconds. This can be expected since the algorithm is based on the 
assumption that there are only a few modes present (Schwartz, 1978), while at increasing 
frequency more modes will appear in the mode separation diagram. A selection of the 
accompanying mode separation diagrams is shown in figure 2.3.8. 
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Figure 2.3.7. Mean and standard deviation of v, for 7 higher modes as a function of 
period. For sake of clarity an offset is introduced along the period axis. Groups of three 
measurements belong to the period of the center of this group. 


The effect of the clean algorithm on the mode separation diagram is that one extra mode 
(mode 4) is resolved, Unfortunately, this has as a consequence that the first and second 
mode give biased results. The PDFT estimator clearly enhances the energy around the 
modes, but if ¢ becomes too small (case d.), energy is smeared out. For this test we calcu- 
lated 14 mode separation diagrams. Each diagram contains 7 higher modes and thus 
represent many different situations of possible mode-mode interference effects and relative 
excitations. This test is therefore thought to be representative for the general situation we 
encounter in practice and not only for the used synthetic data set. The accuracy in the meas- 
urements is not significantly improved when using the PDFT estimator, but the increase in 
the signal to noise ratio is very helpful in identification of the modes. 

For the measurements, presented and discussed in the next chapter, we use a DFT stack first 
to get an impression of the quality of the data and then a PDFT stack with different values 
for €. 
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Figure 2.3.8. Mode separation diagrams at a period of 22.85 seconds for: a. DFT 
stack;b. DFT+clean; c. PDFT stack with €=.1%, box width Ak, , =5% and perturbation 


Ak, ,=1% and d: €=.01%, box width Ak, ,=5% and Ak, , =1%. 


3. ANALYSIS OF NARS DATA 


In this chapter we will present the NARS data that are used for higher mode surface wave 
analysis. Results are presented in the form of a phase velocity dataset for Rayleigh and 
Love fundamental and higher modes. In addition, a shear wave travel time dataset, compiled 
from ISC bulletins, is presented. It is used to constrain the shear velocity in the upper part 
of the mantle. 


3.1 Data 


Figure 3.1.1. Location of earthquakes used in this analysis. 


For optimal analysis of higher modes, earthquakes must be located near the great circle 
through the NARS stations. During 1983,1984 and 1985 three Japanese events and one 
Kurile event (see table 3.1) were recorded in a sufficient number of NARS stations (7 or 
more) to carry out a higher mode analysis. One additional Japanese event (26/5/83) has 
been used that was recorded with good quality records in only 6 stations. The reason to 
include this event is that its location coincides with the event of June 21, 1983 and their 
moment tensor solutions are very similar, so an attempt was made to process these two 
events together in one stack. Surprisingly one event in NW Africa occured on the great cir- 
cle at the end of 1983 and this event was recorded in 7 stations. 

Higher mode surface wave excitation is most effective for a source at intermediate depth. 
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Unfortunately none such event with sufficient magnitude occured in the first three years of 
NARS’s existence. 

Figure 3.1.1 shows the location of the events. The NW Africa (22/12/83) and Kurile 
(24/03/84) events are situated precisely on the great circle through the NARS array. The 
other Japanese events deviate in azimuth with a maximum of 12 degrees for the March 6, 
1984 event. Since this is considerably less than the maximum deviation in azimuth in previ- 
ous studies, which gave good results, we assume that we can neglect the variation of the ini- 
tial phase. This simplifies the application of the stacking technique. The vertical com- 
ponent sections of the analysed events are displayed in figures 3.1.3 to 3.1.5. The horizontal 
components, although less stable, were used to collect a Love wave dataset. 

The data were checked for possible timing errors and polarity reversals by comparison of 
P,S onsets and surface wave polarities along the array. Events recorded during 1983 do not 
contain minute marks on the EW component. In this case the onset times of P and S arrivals 
were calculated using the PEM-C model and compared with the data. Epicenter determina- 
tions have been taken from the ISC bulletin. 

The radial and transverse components of the seismograms of the Honshu/Kurile events are 
displayed in figures 3.1.6 to 3.1.10. The horizontal components of the NW Africa event 
were too unstable to be included in the dataset. 


Table 3.1: List of events used in this analysis. 
The last column gives the number of stations that recorded the event. 


date origin time lon lat depth mb Region # stat. 
26/05/83 2:59:58.8  139.09E 40.48N 16 6.6 Honshu 6 
21/06/83 6:25:26.8  139.10E 41.35N 6 6.5 Honshu 8 
22/12/83 4:11:29.3 13.51W = 11.85N 11 6.3 NW Africa 7 
06/03/84 2:17:21.0 138.92E 29.35N 454 16.1 ~~ &E. Honshu 9 
24/03/84 9:44:02.3 148.12E 44.04N 42 6.1  Kuriles 10 
18/09/84 17:02:40.7  141.64E 34.01N 20 6.5 Honshu 9 


The vertical and radial component of the shallow Japanese/Kurile events show a pulse like 
arrival with a group velocity of approximately 5.2 km/s and dominant period of 40 seconds 
between the predicted SS and SSS (see fig 3.1.2). Nolet (1976) showed that this arrival can 
be caused by a simultaneous constructive interference for both horizontal wavenumbers and 
radial orders, or when group velocities of adjacent modes are well separated but contain an 
independence of frequency in a wide frequency interval. The transverse component shows 
a similar phase at a group velocity of 4.5 km/s. Nakada and Hashizume (1983), in a study of 
the upper mantle beneath the Canadian shield, used similar pulse like signals as higher- 
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Levu | 


mode surface wave data and identified the signal at the transverse component as a S, phase. 
A similar phase on the vertical and radial component with group velocity between 4.6 and 
4.8, called the Rayleigh type S$, phase, can not be found in our dataset. 
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Figure 3.1.2. Vertical, radial and transverse component of event 24/3/84 recorded in 
station NE04. Marked with a star are the pulse like arrivals 
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Figure 3.1.3. Vertical component NARS recordings of the events listed in table 3.1 
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Figure 3.1.4. Vertical component NARS recordings of the events listed in table 3.1 
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Figure 3.1.5. Vertical component NARS recordings of the events listed in table 3.1 
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21/6/83 coordinates: 139.10E 41.35N; depth=6 km;Mb=6.5 
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Figure 3.1.6. Radial and transverse component NARS recordings of the events listed in 
table 3.1 
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Figure 3.1.7. Radial and transverse component NARS recordings of the events listed in 


table 3.1 
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Figure 3.1.8. Radial and transverse component NARS recordings of the events listed in 
table 3.1 
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Figure 3.1.9. Radial and transverse component NARS recordings of the events listed in 
table 3.1 
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Figure 3.1.10. Radial and transverse component NARS recordings of the events listed in 
table 3.1 
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3.2 Measurements 
3.2.1. Higher mode surface wave measurements 


Vertical component recordings of the 6 events of table 3.1 have been analysed in 
order to obtain Rayleigh wave phase velocities. Transverse components could only be used 
for 2 Japanese events. For each event mode separation diagrams of phase- versus group 
velocity have been calculated for 23 frequencies between 20 and 70 mHz. These measure- 
ments are obtained in two runs of 14 frequencies, each run with a different group velocity 
window width, with 5 overlapping frequencies between 30 and 50 mHz. In this way the 
influence of the width of the group velocity window could be checked. Modes, showing up 
as energy maxima, are automatically picked from the diagrams in the neighbourhood of 
theoretical values predicted by the M7 model (Nolet, 1977). Care has been taken that the 
region searched for a maximum is sufficiently large to rule out any bias by the M7 model. 
All diagrams are visually inspected to check whether a real maximum is present and not 
biased by interference effects of neighbouring modes. A first analysis is done without any 
filtering or weighting to get an impression of the data quality. Then several runs with the 
PDFT estimator applying different values for the damping parameter are carried out to 
optimize these results. These data were removed whenever repeated measurements with 
varying damping parameter produced inconsistent results. Application of the stacking tech- 
nique assumes a laterally homogeneous sub-surface structure. The path under study does 
contain some lateral variations (Spakman, 1986; Panza et al, 1980; Nolet et al, 1986d). We 
expect this phenomenon to express itself in a scatter of the measurements, 


Rayleigh waves: 


In order to get an impression of the reliability of the higher mode Rayleigh wave 
phase velocity measurements we calculated vertical component synthetic seismogram sec- 
tions for all 6 events by means of mode summation, using mode 1 to 15 in the frequency 
range of 20 to 80 mHz. For the earth model we used the outcome of a preliminary analysis 
(WEPL1; Dost, 1986), modified by addition of traveltime information in the inversion pro- 
cedure. At depths greater than 600 km the PREM model has been adopted. The synthetic 
sections show the observed pulse like signal at a group velocity of approximately 5.2 km/s. 
Figure 3.2.1 shows a flat part in the curves for mode 7 and 8 at periods between 30 and 40 
seconds at the appropriate group velocity of 5.2 km/s. This confirms the secondly proposed 
nature of the signal in chapter 3.1. 
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Figure 3.2.1, Group velocity as a function of period for the modified WEPL1 model 
used in calculation of synthetic seismogram sections. 


The results of the analysis of the synthetic sections are displayed in figure 3.2.2. These 
measurements have been made without any weighting or deconvolution. From this experi- 
ment we can conclude that: 


[a] 


[b] 


Ic] 


The first higher mode is badly measured, except for the Kurile event (24/3/84). A 
plot of the relative excitations of the different modes, as shown in figure 3.2.3, indi- 
cates that this is due to a bad excitation. Side lobe effects from neighbouring modes 
can strongly influence the measurements. Side lobes coming from the interfering 
modes 5-8 can also generate features in the neighbourhood of mode 1, and can be 
recognized by their high group velocity values. Identification of the modes becomes 
more difficult at higher frequencies where group velocity curves of the various modes 
approach each other. 

The second and third mode can be accurately determined. It is however necessary to 
average phase velocities from all events in order to obtain an unbiased result. 

From mode 4 onward only the high frequency part of the measurements (<25s) gives 
unbiased results. Since we can model the previously described interference pattern 
that causes the bias in the phase velocities, we can try to discriminate between dif- 
ferent earth models obtained after inversion by comparing the waveform of the pulse 
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like arrival with synthetics. 
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Figure 3.2.2. Phase velocity measurements on synthetic seismogram sections The 
different mode branches are labeled. 


This analysis gives us some insight in the reliability of our measurements. We will 
present and discuss the measured phase velocity data set with the results from our synthetic 
experiment in mind. Rayleigh wave measurements that were judged to be reliable are 
presented in figure 3.2.4 and compiled in appendix D. 

The fundamental mode Rayleigh wave is measured up to 50 mHz (20 seconds). The reason 
to exclude higher frequencies is that the penetration depth of surface waves decreases with 
increasing frequency, thus making the surface waves more susceptible to lateral hetero- 
geneity. For the fundamental mode at frequencies higher than 50 mHz this gave rise to an 
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unacceptable scatter in the phase velocity measurements. The same phenomenon is 
observed for the first higher mode. 


EVENT 24/3/84 RELATIVE EXCITATION (MAX MODE O=100) 


20. 30.0 40.0 30.0 
period [s] 


Figure 3.2.3. Relative excitation of mode 1-15 for the synthetics of event 24/3/84 as a 
function of period. Curves for mode 1-9 are marked. 


The fact that the second and third higher mode do not show this large scatter gives us 
confidence that the deeper structure is much smoother. The first higher mode shows a large 
difference in phase velocity between the 21/6/1983 Honshu event and the other events for 
periods between 35 and 50 seconds. This is also the case for mode 2 and 3. As can be seen 
from the data (figure 3.1.3.) the two events in 1983 both suffer from a lack of data in the 
central part of the array. Since these two events have nearly the same location and their 
moment tensor solutions are very similar an attempt was made to process these events 
together in one stack. For the first two modes the phase velocities are reduced and match the 
values for the 1984 events between 30 and 40 seconds. Between 40 and 50 seconds this is 
not the case. This last period range is at the low frequency edge of the interval where we 
can accurately measure higher mode phase velocities with the NARS array. Therefore we 
did not include phase velocities between 40 and 50 second period from the first two higher 
modes into our data set. It should be noted that this is exactly the period range where there 
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are differences in measurements between Cara et al. (1980) and Nolet (1977). Moreover, 
even the values are comparable (see table 3.2). 


Table 3.2 Comparison of phase velocity measurements. 


T [sl vy [km/s] for the 2” higher mode 
Cara et al. (1980) Nolet (1977) 21/6/83 24/3/84 


18/9/84 


For mode 4 to 6 the measurements for periods <25s are thought to be reliable. Between 25 
and 50 seconds strong interference effects play a major role as seen at the synthetics. 
The difference between the preliminary data set (Dost, 1986) and our final dataset is: 


[a] The preliminary data were obtained with a simple version of the PDFT estimator. In 
this version we assumed that all modes are confined to one finite wavenumber band. 
There was no specific information involved about the modes separately. 


[b] Higher mode measurements were restricted to the frequency range of 30 to 70 mHz. 
The new dataset is extended to the lower frequency part (up to 20 mHz for the third 
higher mode). Higher mode phase velocity measurements for frequencies <20 mHz 
resulted in a large scatter because of the limited length of the array. 


[c] Fundamental mode measurements from Seidl (1971) have been adopted in the prelim- 
inary dataset. At present they have been substituted by newly measured values. 
Measurements from 6 to 50 mHz (20-160 s) show a good agreement with Seidl’s 
values. There is however considerable scatter between 14 and 20 mHz (S50 and 70 s). 


[d] In the final dataset S traveltimes are included. This is described in the chapter 3.2.2. 


Love waves: 


An attempt has been made to measure higher mode Love wave phase velocities. Reliable 
transverse components, measured in at least 7 stations, could be obtained for the 24/3/84 
Kurile event and the 18/9/84 Honshu event (see figures 3.1.9 and 3.1.10). The transverse 
components of the 21/6/83 Honshu event were not used, because of their bad coverage (see 
figure 3.1.7). Results for the fundamental and first four higher mode phase velocity meas- 
urements are shown in figure 3.2.7 and compiled in appendix D. For the fundamental mode 
the group velocity is also shown. The measurements are compared with a Love wave phase 
velocity data set by Cara et al. (1980). This data set is based on measurements of one 
Japanese event recorded in western Europe including scandinavia. In addition we compare 
our results with Love wave measurements for North America, based on recordings of two 
Tonga events (Leveque and Cara, 1983). 
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The fundamental mode shows a continuous dispersion curve for the 24/3/84 Kurile event, 
but a jump in phase and group velocity at 25 mHz for the 18/9/84 Honshu event. The cause 
of this jump can be recognized from the data (see figure 3.2.5), bandpass filtered between 20 
and 30 mHz, The Kurile event shows a undisturbed dispersed wavetrain, while the Honshu 
event clearly shows a beating effect at a group velocity of approximately 3.6 km/s, which is 
the measured group velocity coinciding with the jump in phase velocity. This effect can be 
caused by multipathing due to lateral heterogeneity. Comparison with fundamental mode 
phase velocities from Cara et al. (1980) shows that their measurements coincide with the 
results from the Honshu event up to 25 mHz. Unfortunately, their measurements are not 
extended to higher frequencies and we do not know whether they may be influenced by such 
a beating phenomenon. North American phase velocities are consistently lower, except in 
the vicinity of the jump. Again no comparison can be made for frequencies higher than 25 
mHz, since no observations has been made in that frequency region so far. 

The first higher mode shows a large scatter in phase velocity between 25 and 30 mHz. As 
can be seen from the bandpass filtered data (see figure 3.2.6) the beginning of a Sa phase 
can be recognized at a group velocity of 4.45 km/s for the Kurile event, but not for the 
Honshu event. This can be caused by the greater depth of the Kurile event, which is favour- 
able for the excitation of this phase. However, an additional long period feature in the sta- 
tion with the smallest epicentral distance to the Kurile event, just preceding the Sa phase is 
presumably due to an instability of the horizontal seismometers. This will effect our meas- 
urements, but phase velocity measurements without this registration do not produce accurate 
results in this frequency range either. The second and third mode show a continuous disper- 
sion curve. Comparison with Cara et al. (1980) show that their second higher mode phase 
velocities are influenced by interference effects of the second and third mode. The phase 
velocities for Northern America are very close to our measurements. 

Because of the very limited dataset and the anomalous features in the determination of the 
phase velocities just described, we decided not to invert simultaneously for Love and Ray- 
leigh wave phase velocities, but await more accurate measurements. 


Our dataset for the inversion thus consists of: 

[1] Fundamental mode Rayleigh wave phase velocities from 6 to SO mHz. 

[2] First higher mode Rayleigh wave phase velocities from 25 to 45 mHz, 

[3] Second higher mode Rayleigh wave phase velocities from 25 to 70 mHz 

[4] Third higher mode Rayleigh wave phase velocities from 22 to 70 mHz 

[5] Fourth to sixth higher mode Rayleigh wave phase velocities from 45 to 70 mHz. 


Additionally a shear velocity (Sn) data set is used in the inversion. The motivation for and 
compilation of this data set is the subject of the next section, 
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Figure 3.2.4. Fundamental and higher mode Rayleigh wave phase velocity 


measurements. The dashed line indicates average values used in the inversion. 
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Figure 3.2.5. Comparison of transverse components of NARS registrations of event 
24/3/84 and 18/9/84 near the fundamental mode Love wave. The time axis is reduced by 
a factor A/3.45. 
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Figure 3.2.6. Comparison of transverse components of NARS registrations of event 
24/3/84 and 18/9/84 near the first higher mode Love wave. The time axis is reduced by a 
factor A/4.45. Marked by an asterix is an instability in the transverse component 
recording. 
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Figure 3.2.7. Fundamental and higher mode Love wave phase velocity measurements. 
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3.2.2. Traveltimes 


An additional independent source of information about the upper mantle shear velo- 
city structure that can be included in our inversion procedure are Sn traveltimes from 
selected European events. Several interpretations have been given to this Sn wave, e.g. a 
head wave refracted at the base of the Moho, but Mantovani et al. (1977) show that an 
interpretation as a wave that is guided by the lithosphere and travels with a constant surface 
velocity is more consistent with the character of the phase. For shallow depths of foci and a 
model containing a low velocity zone, it was shown that the Sn phase can be interpreted as a 
lid-wave. This means that its energy is confined to the uppermost few tens of kilometers of 
the mantle. 

Concerning the European continent there have been previous studies that incorporated shear 
wave traveltimes. 

Lehmann (1961) compiled S readings and was able to demonstrate the existence of a low 
velocity layer in Europe, followed by a rise at 220 km depth. Her dataset consisted of ISS 
traveltimes and personal readings from four shallow events (two Italian, one Azores and one 
Greek) and 6 deep Rumanian earthquakes (130 km depth). 

A difference in structure between western and eastern Europe was assumed in the study by 
Mayer-Rosa and Mueller (1973) who examined P and S times together with fundamental 
mode Rayleigh waves. Their conclusion was that in the upper 200 km a clear distinction 
could be made between the two regions. 

Regional S arrivals are difficult to read and may show considerable scatter. Especially too 
large arrival times are to be expected due to a larger amplitude SSn phase that arrives after 
the Sn phase (Kennett, 1985). For selected events and stations that sample the structure 
under NARS we will try to infer a traveltime curve using ISC bulletins as a database. In the 
present compilation we will take care to minimize the influence of lateral heterogeneity and 
therefore neglect events in eastern Europe as well as events in Algeria and southern Spain 
(Spakman, 1986; see also Panza et al (1980) for anomalous Sn velocities in this last region). 
The remaining events and reduced traveltime curve are shown in figure 3.2.8. 
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24/10/76 
6/4/77 
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22/6/82 42.89N/1.94W 10 = 19:50:21.5 
29/7/82 60.31N/1.93E 10 = 00:17:00.2 

8/3/83 59.67N/S.20E 33 = 18:43:54.0 

4/8/83 50.38N/4.33E 10 = 07:08:24.4 

11 9/8/83 50.40N/4.29E 10 = 01:32:34.8 
12 8/11/83 50.73N/5.31E 5 00:49:31.5 
13 8/11/83 50.60N/5.38E 8 = 02:13:21.2 
14 3/4/84 SO.S8N/249W 12 = 14:27:11.2 
15 15/4/84 52.22N/.31W = 25 21:05:34.4 


T -—20A [s] 


Figure 3.2.8. Top:list of used events and great circle paths. Bottom: Reduced Sn traveltime as a 
function of epicentral distance A. Solid line represents the inferred traveltime curve, dashed lines 
the assumed standard deviation. Data from two french events, event 6 and 7, are marked with a 
cross in stead of a triangle and are not used in construction of the traveltime curve. 
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The dataset contains mainly events in the North sea, United Kingdom, Belgium, the Nether- 
lands and stations in France. The inferred traveltime curve is compatible with the majority 
of the data. A scatter in the measurements between 5 and 12 degrees is mainly due to 2 
French events, number 6 and 7. In general regional French events give a large scatter in Sn 
traveltimes and are thus omitted. An average standard deviation in the traveltime curve of 
+2.5 s, as indicated in figure 3.2.8, is thought to be representative for our dataset. 

One of the events, number 16, was also recorded by the NARS array in stations NEO2 to 
NEO7. This enabled us to check wether the Sn onset is clearly visible on a broad-band 
instrument or not. Vertical component registrations are shown in figure 3.2.9 
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Figure 3.2.9. NARS registration of the 19-7-1984 United Kingdom event. 


S phases are clearly observable at epicentral distances greater than 700 km. For shorter dis- 
tance the onset is difficult to pick. These measurements have been incorporated in our com- 
pilation of regional traveltimes and are consistent with ISC determinations. 
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4. INVERSION 


The dataset is inverted using the generalised inversion technique (Wiggins, 1972) 
with a sharp eigenvector cut off. In the inversion a trade-off exists between depth resolution 
and variance in the model parameters. The number of eigenvectors used in the inversion 
determines this trade-off. S velocity, v,, and density, p, were inverted simultaneously after 
weighting them with a priori uncertainty. P velocity, v, , is coupled to v, with a fixed Pois- 
sons ratio. This is convenient since only the fundamental mode is weakly sensitive to v,,. 
One should keep in mind that surface waves are only sensitive to changes in the gradient of 
the density, not to absolute changes (Nolet, 1977). As an extra constraint, no change in sur- 
face gravity is allowed to result from the inversion. As a measure of the fit of a model to the 
data, relative residuals are calculated. This is the average root mean square of the differ- 
ence between observed and calculated values (6) for the data divided by the standard error. 
Thus, for N data & corresponds to (x7/N)"? with x? the more familiar statistical quantity 
Since N >>K , the number of eigenvectors used in the construction of the model, the number 
of degrees of freedom is approximately N. Thus, when €<1 a model has been chosen that 
fits the errors in the data rather than the trend in the data. On the other hand, if £>1 a poor 
fit between the data and the model is obtained. Therefore we should aim at a relative resi- 
dual E=1. Of course we might introduce at this stage sophisticated statistical techniques to 
decide whether the addition of an extra eigenvector is warranted by the data. We have, how- 
ever, found it much more convenient to apply subjective judgement to decide whether we 
are near enough to €=1 and to avoid model features of doubtful reality. We shall make up 
for this by applying resolution analysis a posteriori. Resolution is calculated by introducing 
constant velocity changes over thick layers, thereby forcing the system to be overdeter- 
mined, The inversion then gives the variance in the model averaged over the specified layer 
thickness. 

In the inversion procedure a starting model is required. One should be cautious not 
to influence results with a priori assumptions about e.g. a low density- or velocity channel in 
this starting model. The M7 model is an isotropic model. Cara et al. (1980) observed a pos- 
sible anisotropy under western Europe. They based their observation on the incompatability 
of Love and Rayleigh wave phase velocities in this region. Mitchell (1984) mentions as 
possible alternative explanations for a Love/Rayleigh incompatability a laterally varying 
structure or a failure to include a sufficient number of layers in the inversion process. Since 
there is at the moment no proper Love wave dataset at hand for the west-European platform, 
we will interprete our data assuming an isotropic model. 

Nolet’s data sample the Scandinavian shield and the west-European platform, while 
our new data sample only the platform structure. Inversion of the two datasets, using the 
same starting model enables us to construct a shield model. Differences in shear velocity 
and density can thus be investigated. 
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4.1 Inversion for the west-European platform 


In a preliminary inversion (Dost, 1986), using data obtained by an estimator that 
only takes into account the fact that all surface wave modes occur in one limited 
wavenumber band, we took a modified M7 model as the starting model. The M7 model 
does not exhibit sharp discontinuities near 400 and 670 km depth, since the higher mode 
dataset used in constructing this model could not distinguish between a sharp or smooth 
feature at those depths. The shear velocity under the Moho down to 220 km depth was set 
to a constant 4.45 km/sec and density to 3.40 gr/ccm. We restricted density variations in this 
depth interval to 0.2 gr/ccm in order to obtain realistic results. In the preferred model 
(WEPL1) the most remarkable features are a broadening of the low velocity zone and the 
disappearance of a Lehmann discontinuity. In addition a high shear velocity lid between 80 
and 120 km was found to correlate with a density high. In the present inversion with a more 
refined dataset we will investigate the influence of the starting model on the inversion 
results. First the new dataset is inverted with the same smooth starting model as WEPL], 
then we will investigate the influence of sharp discontinuities at varying depths on the inver- 
sion result. An overview of models used in this chapter is given in table 4.8. 


Inversion with a smooth starting model 


The result of our inversion using the modified M7 model as starting model (WEPL 1a) is 
shown in figure 4.1.1. In the same figure WEPL1 is reproduced. Incorporation of travel 
times results in a more pronounced high shear velocity lid and low shear velocity layer. The 
Lehmann discontinuity at 220 km depth is absent in both models and average shear velocity 
between 350 and 650 km is in WEPL1a 0.1 km/s higher than WEPL1. Even if we assume 
the Lehmann discontinuity to be present in the starting model it is removed in the inversion 
procedure. Regarding the density it is remarkable that the density high at 100 km depth 
disappears completely and the density high around 200 km depth remains. The question 
rises whether this density feature is resolved or not. To answer this, resolution is calculated 
at specific depths and given in table 4.1 and 4.2. The standard deviation of the shear velo- 
city v, is comparable for both models. The standard deviation of the density shows a totally 
different behaviour. The WEPL1a model has a resolution far worse than WEPL1. 
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Table 4.1: v, resolution 
WEPLIa WEPLI1 
depth Vs stdev layer thickness | depth Vy stdev layer thickness 
[km] [km/s] _ [km/s] [km/s] __ [km/s] 


Table 4.2: p resolution 
WEPLIia WEPLI1 

depth p stdev _ layer thickness | depth p stdev layer thickness 

[km] [gr/ecm] _[gr/ccm] ({km] [gr/ccm] [gr/ccm] [km] 


100 =. 3.358 .160 80 
180 3.520 .220 80 
260 3.455 


440 3.796 


density [gr /ccm];V, [km/s] ; 


200. depth km] 400-0 600.0 800.0 
Figure 4.1.1. Density (lower curve) and shear velocity as a function of depth. 


This means that a density high at approximately 200 km depth is not resolved and as a 
consequence € will not change when we smooth the density curve. Calculation of the model 
fit to the data for a model with a smooth density curve between 100 and 300 km depth and 
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WEPLIa, however, points towards the existence of the density high: & is reduced from a 
value of 1.69 to 1.38, which is significant. This disagreement between resolution calcula- 
tions and the model fit to the data can be caused by the fact that the linearity assumption in 
the inversion is violated. Experiments show that the thickness of the high velocity lid deter- 
mines the resolution in the density: a smaller lid thickness results in an improvement in 
resolution. It is interesting to note that Lerner-Lam (1985) reports that lid thickness can be 
substantially reduced by allowing a transverse isotropic parameterization of velocity in the 
inversion. In our present inversion this is not accounted for. 

Finally, an experiment has been carried out to test the correlation between y, and p at 100 
km depth, as suggested by WEPL1. Therefore we coupled v,, v, and p by Poisson’s ratio 
and Birch’s law in the inversion. As a result the shear velocity exhibits again the high velo- 
city lid, but density does not show a corresponding high as in WEPL1. 

The model fit to the data shows in all cases a trade-off between the fundamental mode data 
from 90 to 160 s and higher mode data. Mode 5 can not be adequately fitted, the measured 
phase velocities are always too low. This bias can be due to interference effects with the 4” 
and 6” higher mode. 


Inversion with a starting model containing sharp discontinuities 


For a first inversion with a model containing sharp discontinuities we took the modified M7 
model down to 360 km depth and added for depths >360 km the PREM model (Dziewonski 
and Anderson, 1981). In a first experiment the influence of a possible Lehmann discon- 
tinuity is investigated. Three starting models are used in the inversion: one with a smooth 
gradient between 200 and 300 km depth, one with a strong gradient between 220 and 240 
km depth and one with a first order discontinuity at 220 km. Results are displayed in figure 
4.1.2. All models, build with 9 eigenvectors, fit the data equally well (6 =1.44, 1.44 and 
1.42). These values are comparable with results from inversions with a smooth starting 
model. We can conclude that the existence of a Lehmann discontinuity is not resolved with 
our current data set. In general, the pronounced density high at 200 km depth, as present in 
inversions with a smooth starting model, disappears. An increase in density from 40 to 100 
km depth is resolved for all models. For a starting model with a constant density from 40 to 
300 km depth, the density gradient is always introduced and results in a significant improve- 
ment of the fit to the data. This agrees with our WEPL1 model, although the density gra- 
dient is smoother now. Variations in thickness of the high velocity lid do not affect resolu- 
tion in density as was the case in the inversion with a smooth model. Models with a Leh- 
mann discontinuity exhibit a thicker high velocity lid and a more pronounced low velocity 
layer. These models have a better fit to the traveltimes than less pronounced models and are 
thus preferred. The ’400’ km discontinuity starts for all models at higher shear velocity 
values than PREM. The model without a Lehmann discontinuity exhibits a jump in shear 
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velocity that is about 1.8 times larger (Av,= 0.3 km/s) than PREM (Ay,= 0.17 km/s) and 
approximately 1.3 times larger than the two models containing a Lehmann discontinuity 
(Av, = 0.23 km/s). To our knowledge, the only published global or regional model with a 
comparable velocity contrast at the 400’ km discontinuity is the EPR-B model for the man- 
tle under the East-Pacific Rise (Wielandt and Knopoff, 1982). However, the depth of the 
discontinuity in this model is much larger (465 km). The average shear velocity of all three 
models between 400 and 600 km is significantly higher than PREM. Resolution of v, and p 
is approximately the same as in WEPL1a , see table 4.3 and 4.4. 


Table 4.3: v, resolution 
The second column (WEPL1.Ieh) gives values for a model 
with a strong gradient as Lehmann discontinuity. 


WEPL ia WEPL1.leh 
depth Vs stdev layer thickness | depth V, stdev layer thickness 
{km/s} [km/s] [km/s] 
4.605 
4.128 
: 4494, 
340 4.762 .083 80 340 4.717 .078 80 
440 4.954 072 120 440 5.158 — .089 120 


Table 4.4: p resolution 
The second column (WEPL1.leh) gives values for a model 
with a strong gradient as Lehmann discontinuity. 
WEPLIia WEPLI.leh 
depth p stdev layer thickness | depth p stdev _ layer thickness 
[km] [gr/ccm] _[gr/ccm] km] [gr/ccm] _[gr/ccm] 


The disappearance of the density anomaly at 200 km depth and the disappearance of the 
dependence of density resolution on lid thickness stress the importance of the existence of 
discontinuities in the starting model. No explicit assumptions about anisotropy or a density 
high at 200 km depth are needed to obtain an acceptable model. The departure from smooth 
models in the linearized inversion may introduce non-linear effects. The presented resolu- 
tion calculations should therefore be interpretated as a rough order of magnitude calculation. 
In conclusion we can say that, although the Lehmann discontinuity is not resolved, a smooth 
gradient between 200 and 400 km depth leads to a Jarge jump in shear velocity at the ’400’ 
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km discontinuity compared to a reference earth model like PREM and more regional 
models. 


OL | NO LEHMANN DISCONTINUITY 4 


mf STRONG GARDIENT AS LEHMANN DISCONTINUITY 
re arte FIRST ORDER LEHMANN DISCONTINUITY 
= PREM 


density [gr/ccm] ; V, [km/s] 


! 2 | 1 = t 
200.0 400.0 600.0 800.0 


depth [km] 


Figure 4.1.2. Density (lower curve) and shear velocity as a function of depth for a 
starting model with discontinuities at 400 and 650 km depth. Inversion results for three 
different starting models are compared with PREM. 


Now we can investigate the influence of the depth of the two major discontinuities in the 
inversion. The depth of the 400 km discontinuity is best determined by body-wave studies 
and ranges from 380 to 420 km depth. In general the sharpness of this discontinuity is not 
well determined, except in the case of northern Australia where a first order discontinuity at 
406 km has been found (Leven, 1985). We investigate three models with a first order 
discontinuity at 380, 400 and 420 km depth respectively as well as a model with a strong 
gradient between 380 and 420 km depth. Results for the first three models are displayed in 
figure 4.1.3. The discontinuity located at 380 km depth has the effect that the Lehmann 
discontinuity disappears and that the jump in v, at the discontinuity becomes much larger 
(.28 km/s) than in PREM. The other two models are very similar to each other, except in 
the immediate vicinity of the discontinuity itself. The jump in shear velocity at the discon- 
tinuity is in both models approximately .22 km/s. This value is comparable with existing 
regional models like SNA and TNA (Grand and Helmberger, 1984). The gradients in shear 
velocity between the ’400’ km discontinuity and 600 km depth are approximately equal for 
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the three models. There is no resolved difference in density gradient for the three models. 
Experiments in the form of an inversion with a strong gradient in shear velocity and density 
between 380 and 420 km depth gave similar results as the models with a discontinuity at 
400 and 420 km depth. This observation shows again the inability of the surface wave data 
to discriminate between a strong gradient and a first order discontinuity. The influence of 
the location of the 650 km discontinuity on more superficial gradients in shear velocity or 
density of the model is negligible. 

There is a difference in width of the high velocity lid and absolute depth of the low velocity 
layer between inversion results in figure 4.1.2 and 4.1.3. As stated before it was found that 
models with a greater thickness of the high velocity lid and more pronounced low velocity 
layer give the best fit to the traveltime data and are thus to be preferred. The use of models 
with a more pronounced velocity contrast and thicker high velocity lid in stead of the 
models shown in figure 4.1.3. has no influence on our present conclusions. 
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Figure 4.1.3. Density (lower curve) and shear velocity as a function of depth for a 
starting model with the ’400 km’ discontinuity at 380, 400 and 420 km depth. Inversion 
results for three different starting models are compared with PREM. 


For a definite choice of an average platform model we include constraints from body wave 
studies. Paulssen (1987 and pers. comm.) studied P and S waveforms from European events 
using NARS recordings and WKBJ modelling. She found that waveforms from different 
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source-receiver paths could not be fitted with one homogeneous model down to a depth of at 
least 450 km and only bounds on the depth of the ’400’ km discontinuity can be roughly 
specified. A minimum depth of the ’400’ km discontinuity of 400-420 km is constrained 
from this modelling. 

Our final model for the west-European platform (WEPL2) thus consists of a high shear 
velocity lid from 80 to 140 km depth, a pronounced low velocity channel from 140 to 240 
km depth and a smooth gradient in velocity between 240 and 400 km depth. The two major 
discontinuities are placed at 400 and 650 km depth. Density shows a steep gradient between 
60 and 140 km depth, followed by an approximately constant value down to 300 km depth. 
At greater depth density follows PREM closely. This model is shown in figure 4.1.4 and 
listed in table 4.5. Resolution is indicated in table 4.6. In all inversions presented in this 
chapter there is a trade-off between the misfit of fundamental mode data at 90 to 160 s and 
higher mode data as noted already in the inversion with a smooth starting model. A period 
of 90 s for the fundamental mode means a penetration depth of approximately 180 km, the 
deepest part of the low velocity zone. The region that is of major importance in this incom- 
patibility between fundamental and higher modes is thus the region around the possible 


Lehmann discontinuity which is badly resolved in the inversion. The bias of the 5“ mode 
remains also present. 


Table 4.5: WEPL2 


depth v 


[km] [km/s] 
0-5 3.100 
5-20 3.509 
20-29 3.829 
29 4.263 
40 4.323 
60 4.382 
80 4.543 
100 4.622 
120 4.622 
140 4.520 
160 4.178 
180 4.117 
200 4.130 
220 4.244 
240 4.414 
260 4.453 
280 4.505 
4.557 
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Figure 4.1.4. Density (lower curve) and shear velocity as a function of depth for 
WEPL2. For comparison M7 and PREM are shown. 


Table 4.6: resolution WEPL2 
depth Vs stdev layer thickness | depth re) stdev _ layer thickness 


{km/s] _ [km/s] [km] [gr/ccm] _ (gr/ccm] 


4.2 Inversion for the Scandinavian shield 


We can compare our new dataset sampling the west-European platform with Nolet’s 
dataset for western Europe including Scandinavia in order to detect differences in structure 
between the two regions. Although our measurements have been obtained for much higher 
frequencies, this comparison can give information on gross differences. Therefore we com- 
pared smooth models inverted with only a limited amount of eigenvectors. Nolet’s phase 
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velocity data were re-inverted using the same starting model as in the WEPL2 inversion. A 
smooth model could be obtained with 7 eigenvectors and shows basically the same features 
as M7. A minimum in density around 220 km depth is retained. The new data were re- 
inverted to obtain a smooth 7 eigenvector model for comparison. Traveltimes have been 
excluded. The shield model has been calculated by means of the expression: 


Mt, A, = (Ap +A, heap — Ap itty (4.1) 


where A, and A, stand for the path length over the shield and the platform structure respec- 
tively and 7, 77, and 77,,, denote shield, platform and averaged model values. The result- 
ing model is called SCSH2 and is tabulated in table 4.7. The same procedure has been 
applied to the preliminary dataset by Dost (1986), and resulted in a difference model 
SCSH1. A comparison between SCSH1, SCSH2 and WEPL2 is shown in figure 4.1.5. A 
constant shear velocity value down to 300 km depth that was found in SCSH1, is still 
present in SCSH2. The latter model shows more fluctuations between 220 and 300 km 
depth, but this can be attributed to the region where resolution becomes bad in WEPL2 and 
the possible existence of a Lehmann discontinuity is not resolved. This observation that the 
low velocity channel is absent under the Scandinavian shield is also supported by findings 
of Calcagnile (1982) from regionalisation of fundamental mode Rayleigh wave phase velo- 
cities. Between 400 and 500 km depth a difference in shear velocity between WEPL2 and 
both SCSH1 and SCSH72 is observed. Concerning the density, a difference between 200 and 
320 km depth is observed. This density difference is on average 0.20 gr/ccm. As we noted 
before resolution calculations can only be used as a rough estimate. In these model differ- 
ence calculations this becomes even more troublesome. However, since the density differ- 
ence is spread out over at least 120 km depth and the estimated resolution for WEPL2 at 
260 km is 0.21 gr/ecm over 80 km depth, it is acceptable to conclude that there is a strong 
indication for a difference in density between the Scandinavian shield and the west- 
European platform between 200 and 350 km depth. The interpretation of the models 
presented in this chapter in terms of mineralogy is the topic of the next chapter. 
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Figure 4.1.5. Density (lower curve) and shear velocity as a function of depth for 
SCSH2. For comparison SCSH1 and WEPL2 are shown. 


Table 4.7: SCSH2 


depth Vy p depth Vs 


p 
[gr/ccm] [km] [gr/ccm] 
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Table 4.8 Overview of models mentioned in the text. 


model data used # eigenvectors remarks 

M7 Fundamental (6-35 mHz) and higher ? Separate inversion 
mode (mode 1-6;10-35 mHz) Rayleigh for shear velocity 
wave phase velocities for western and density. 
Europe and Scandinavia (Nolet, 1977). 

WEPL1 Fundamental (6-50 mHz) and higher 10 Simultaneous 
mode (mode 1-6;30-70mHz) Rayleigh inversion for shear 
wave phase velocities for the west- velocity and density 
European platform (Dost, 1986). 

IWEPL2 Fundamental (6-50 mHz) and higher 9 Simultaneous 
mode (mode 1-6;22-70mHz) Rayleigh inversion for shear 
wave phase velocities for the west- velocity, density and 
European platform. In addition a shear shear wave 
wave traveltime dataset is compiled. traveltime data. 


Discontinuities are introduced in the 

starting model (This study). 

SCSHI1 Constructed as a difference model from 8 see WEPL1 
an inversion of the dataset used in 

WEPLI and the M7 dataset using 

identical starting models (Dost, 1986). 

SCSH2 Same as in SCSH1, only WEPL! is 7 see WEPL2 
replaced by WEPL2 and traveltimes are 

excluded (This study). 


WEPE la Inversion of the dataset of WEPL2 with 10 none 
a smooth starting model (This study). 
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5. INTERPRETATION AND CONCLUSIONS 


For an interpretation of the models derived in the last chapter in terms of mantle 
mineralogy we need information about elasticity parameters and their pressure and tempera- 
ture derivatives. Recent progress in the measurement of these parameters in high pressure 
laboratories now make it possible to test hypotheses about the upper mantle structure. 
These hypotheses are mainly based on arguments related to the average mineralogy (pyrol- 
ite and piclogite models), but may include arguments concerning mechanical properties 
(tectosphere model or lithospheric doubling concept). We will give a brief review of the 
major hypotheses, compare our shear velocity and density models with models constructed 
from recently published laboratory data and discuss the implications for the west-European 
platform and Scandinavian shield structure. An interpretation of preliminary results 
(WEPL1 and SCSH1; Dost, 1986) in terms of mineralogy was published in Nolet et al. 
(1986b). 


5.1 Mineralogical models 


Table 5.1 gives an overview of the principal minerals in the upper mantle. Indicated 
are abbreviations that will be used in the text. 


Pyrolite 


Ringwood (1975, 1979, 1982) proposes a model for the upper mantle that is based 
on a homogeneous undepleted mantle mineral assembly: pyrolite. Its composition, as used 
in different tests mentioned in this chapter, is shown in table 3.2. Pyrolite acts as the primi- 
tive source material of modern basalts, that are formed upon partial melting. Different 
basalt types are the result of a difference in the percentage of partial melting (this percen- 
tage is for nephelinites 1-5%, alkali basalts 5-10%, tholeiites 15-25% and komatiites 30- 
60%). The remaining fractions after partial melting are garnet and spinel lherzolite, com- 
plementary to nephelinites and alkali basalts, harzburgite, complementary to tholeiites and 
dunite, complementary to komatiites. The discontinuities at 400 and 650 km depth are inter- 
preted as phase changes. The ’400’ km discontinuity is interpreted as the transition from 
olivine to modified spinel (B-phase) accompanied by a increasing solid solution of pyroxene 
in garnet. The transition zone, the region between 400 and 650 km depth, is mainly con- 
sistent with a modified spinel plus garnet composition with minor clinopyroxene. The latter 
phase exists down to 500 km depth. At 550 km depth the modified spinel transforms to 
spinel (y-phase) and this breaks down to perovskite and periclase between 650 and 700 km 
depth. Garnet breaks down to the ilmenite structure starting at 600 km depth. A subducting 
slab, mainly consisting of a basaltic upper layer underlain by a harzburgite layer, sinks to 
approximately 600 km depth and remains denser than the surrounding mantle. Below 650 
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km depth the former basaltic layer remains denser, but the harzburgite layer becomes rela- 
tively buoyant and material gathers at the 650 km discontinuity. This material is heated by 
thermal conduction and the former basaltic crust is subjected to partial melting. As a result 
the former harzburgite is fertilized, becomes buoyant and starts forming diapirs into the 
upper mantle. The residual former basaltic layer sinks as large blocks into the lower mantle 
in a dense perovskite structure. 


Piclogite 


Anderson (1979a), Bass and Anderson (1984) and Anderson and Bass (1984, 1986) discuss 
the piclogite model in detail. A general discussion is also presented by Oxburgh (1980). 
The piclogite model is based on a chemically inhomogeneous mantle composition. It is 
assumed that an original melt layer or magma ocean was present. Then the scenario is 
envisaged where low density olivine crystallizes and becomes concentrated in the shallow 
mantle in peridotite. Garnet crystallizes later and sinks, while the residual fluid freezes to a 
clinopyroxene-garnet rich mixture with some olivine. This residu is called piclogite (or an 
olivine eclogite) and is in composition close to komatiite. Pyroxene dissolves with increas- 
ing depth in the garnet and at the end of the transition zone (550-650 km) garnetite is 
formed. Eclogite will accumulate at the 650 km discontinuity, as it cannot penetrate it 
because of the high density of A!,04- poor assemblages beneath the discontinuity. In this 
way a Stratification of the upper mantle is achieved. From top to bottom the upper mantle 
thus consists of: basalt, peridotite and eclogite. Basalt injected into the mantle during sub- 
duction inverts at shallow depths to eclogite and sinks rapidly through the mantle. Ander- 
son (1979a,b) interpretes the 220 km discontinuity as the boundary between the peridotite 
and eclogite layer and the 400 km discontinuity as the boundary between eclogite and gar- 
netite. More recently Liu (1980) and Bina and Wood (1984) have argued that the eclogite 
to garnetite transition would not produce the observed 400 km discontinuity. As a result 
Anderson and Bass (1984, 1986) and Bass and Anderson (1984) moved the boundary 
between peridotite and eclogite down from 220 km to 400 km depth. Most remarkable is 
the difference in piclogite composition between the Bass and Anderson (1984) and Ander- 
son and Bass (1984) paper (see table 5.1). A comparison between compositions used by 


different authors is complicated by the use of both weight and volume percentages by these 
authors. 
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Table 5.1 Principal minerals in the upper mantle. 


Olivine (Ol) (B= modified spinel, y=spinel) [Mg ,Fe],SiO, 

Orthopyroxene (Opx) (Mg ,Fe)SiO3 

Clinopyroxene (Cpx) (Mg ,Fe )CaSi20 ¢ 

Pyrope-garnet (Pgt) (Mg Fe )3Al,Si30 12 

Garnet-majorite (Mgt) [(Mg Fe )3Al 2Si30 12)1-, [(Mg Fe )4Si49 121, 
Perovskite (Pv) (Mg ,Fe)Si0 3 

Ilmenite (11) (Meg Fe)SiO 5 

Periclase (Pc) (Mg Fe)O 


Jadeite (Jd), commonly in solid solution NaAlSi,0 ¢ 


with Cpx 


Grossularite (Grs), commonly in solid CazAl2Si30 12 
solution with Pgt 
Garnet lherzolite Pgt+Opx+Cpx+Ol 
Spinel lherzolite  Peridotite Sp+Opx+Cpx+Ol 
Harzburgite O1l+Opx 
Dunite Ol 
Eclogite Gt+Cpx 
Table 5.2 Comparison of models used in hypothesis testing. 
Pyrolite Piclogite 

Comp. | Bass & And. And. & Bass Weidner | Bass& And. And. & Bass Weidner 

[wt %] [wt %] [vol %] [wt %] [wt %] [vol %] 
Ol 57 57 61 16 16 16 
Opx 17 17 15 3 6 3 
Cpx 12 12 10 44 33 45 
Gt 14 14 14 37 45 36 

Hypothesis testing 


In a recent series of articles, partly reviewed by Anderson and Bass (1986), the 
pyrolite and piclogite models are tested against global seismological results (the PREM 
model of Dziewonski and Anderson, 1981). As a data base an internally consistent dataset 
of laboratory measurements on a pyrolite composition by Akaogi and Akimoto (1977, 1979) 
is widely used. A review by Jeanloz and Thompson (1983) gives valuable information on 
existing data and uncertainties. The only laboratory measurements on a piclogitic composi- 
tion is made by Irifune et al. (1986). They studied an alkaline-poor olivine tholeiite. 
Acoustic velocities are calculated at depth by projection along adiabats. This assumption of 
adiabaticity is thought to be reasonable for the deeper part of the upper mantle. In terms of 
thermodynamics the mineral assemblage is heated isobarically up to the so called ’foot’ to 
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the adiabat temperature. From there it is adiabatically compressed to the desired pressure 
(Weidner, 1986). 

Bass and Anderson (1984) and Anderson and Bass (1984) compared in this way a pyrolite 
and piclogite composition. As reference model the PREM model has been used. They con- 
cluded that both the pyrolite and piclogite models are acceptable down to the 400 km 
discontinuity, but that for the transition region between 400 and 650 km depth the pyrolite 
model generates too high shear velocity values, while piclogite is appropriate. Anderson and 
Bass (1984) show that for shields the velocities fall near the 1400°C adiabat for depths 
greater than 200 km. Other models, e.g. tectonic rise, approach this adiabat at much larger 
depths ( near 400 km). In answer to their study Weidner (1985; 1986) comes to the opposite 
conclusion, namely that the pyrolite model is acceptable, and points out that this incon- 
sistency is due to different assumptions concerning phase equilibria for the pyroxene-garnet 
transformation, elastic properties of the majorite phase and pressure and temperature deriva- 
tives of the elastic properties of the high pressure phases. In his study a 1400° C foot to the 
adiabat has been used as well. Irifune et al. (1986) conclude that, based on measured pro- 
perties of a piclogite like composition, the piclogite model can not explain the velocity jump 
at 400 km depth and the high gradient in the transition zone. 

The most important parameters that lead to the present controversy are the jump in shear 
velocity at the 400 km discontinuity, the velocity gradient in the transition region (400-650 
km depth) and the absolute level of velocity between 300 and 400 km depth. In the follow- 
ing we will describe these in detail. 


Absolute velocity 300-400 km depth 


In this region all calculated models give more or less the same results. Here the assumption 
of the 1400°C foot to the adiabat is only valid for shields. 


400 km discontinuity 


Bass and Anderson (1984) assume in their model that at 400 km depth all olivine is 
transformed to modified spinel ( B phase) and that all orthopyroxene is already dissolved in 
garnet as majorite. Their pyrolite model therefore creates a jump in velocity that is twice 
the observed value and only piclogite, containing less olivine, is able to produce the correct 
velocity jump due to the additional pyroxene to garnet transition. Bina and Wood (1984) 
demonstrate, both on experimental and thermodynamical grounds, that the eclogite-garnetite 
transition does not exhibit any discontinuity in bulk sound velocity from 11 to 23 GPa. 


86 


Weidner (1985) argues that the magnitude of the 400 km discontinuity is not only sensitive 
to the amount of olivine present in the mantle, but also to the pressure derivative of the 
shear modulus of the high pressure phase that is not yet measured. The observed velocity 
and density jump is found to be consistent with an olivine content of 40-65 volume %. For 
models with a smaller olivine content an accompanying chemical discontinuity should be 
present. The contradiction with findings by Lees et al. (1983) are attributed to incorrect 
elastic parameters for spinel and modified spinel and differences with Bass and Anderson 
(1984) to phase relations and properties of the garnet phase. In a more elaborate paper on 
the subject Weidner (1986) argues that there are observations that are incompatible with 
Bass and Anderson’s conclusions. One critical assumption that comes under attack is the 
assumption that Ca-poor pyroxene dissolves completely into the garnet phase around 400 
km depth. The only possibility for a complete dissolution of Ca-poor pyroxene is in case 
Akaogi and Akimoto (1977; 1979) overestimated their pressure. Irifune et al. (1986) report 
that this is the case. A direct cause of the overestimation of pressure is that there should be 
a considerable solution of pyroxenes into garnets at depths as shallow as 200 km. The recent 
discovery by Moore and Gumey (1985) of pyroxenes in solid solution with garnets as inclu- 
sions in diamonds from kimberlites supports these conclusions. 

Observations that can not be explained are coexisting garnet, modified spinel and Ca-poor 
pyroxenes at a pressure of 145 kb and spinel with Ca-poor pyroxene at 205 kb (Akaogi and 
Akimoto, 1979). 

The experiments carried out by Irifune et al. (1986) on a piclogitic composition lead to the 
conclusion that not only the Ca-poor-pyroxene dissolves completely into the garnet phase at 
400 km depth, but also the clinopyroxene. This gives us a pyroxene free garnetite composi- 
tion at the 400 km discontinuity, which results in a too large jump in velocity and density for 
Anderson and Bass’s eclogitic and piclogitic compositions. 


Transition region 400-650 km 


Anderson and Bass (1986) assume that this region is dominated by a clinopyroxene-garnet 
composition. At 550 km depth, 60% of the clinopyroxene is dissolved in the garnet phase 
as majorite and coexists with spinel and orthopyroxeen-majorite. At the bottom of the tran- 
sition layer all pyroxene is present as majorite. They claim that only this composition, with 
clinopyroxene gradually changing to majorite can explain the low velocities and high gra- 
dients in the transition region. Weidner (1985, 1986) demonstrates that his pyrolite model 
does account for the observed high gradient and concludes that the pyroxene to garnet and 
excess Ca-poor-pyroxene to spinel plus stishovite transformations give the main 
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boundaries between modified spinel and spinel remain linear, but that the above mentioned 
interpretation of the 550 km discontinuity has to be rejected. The same author proposes a 
phase transition of pyroxene to a garnet-like structure. 


5.2 Models involving mechanical properties 


The models described in the last section are based on an average Earth structure without 
lateral variations. In this chapter we will discuss two models that do account for differences 
between geological provinces, 


Tectosphere model 


Jordan (1975, 1978, 1981) tackled the problem of the contrast in upper mantle properties 
between old continents and old oceans. Application of the thermal-boundary layer model 
(as developed in plate tectonic theory) to continents preclude the existence of differences in 
chemical, mechanical or thermal properties beneath 100-150 km depth. However, from 
seismological observations of ScS and multiple ScS travel times (Sipkin and Jordan, 1975, 
1976, 1980) this was questioned and the tectosphere model presented as an alternative. In 
the tectosphere model a thick ’root’ or chemical boundary layer is present that stabilizes a 
more than 200 km thick thermal boundary layer against convective disruptions. The root 
consists mainly of peridotite, depleted in basaltic components. Depletion can be caused by 
10-20% partial melting. The melt phase migrates upward and as a result the root becomes 
less dense. This effect will not only compensate for the temperature, but also stabilize 
against convection. This hypothesis has been rejected by e.g. Anderson (1979b), because of 
a lack of convincing seismological evidence that continents and oceans differ in structure 
beneath 200 km depth. 


Lithospheric doubling 


Lithospheric doubling is defined by Vlaar (1982, 1983) as horizontal subduction of young 
(<30 Ma) oceanic lithosphere, obducted by a continental or oceanic lithosphere. It is pro- 
posed to be a general phenomenon in recent plate tectonics and was shown to provide a 
succesful scenario for the Alpine orogeny (Vlaar and Cloetingh, 1984). Lithospheric dou- 
bling would imply a strongly stratified upper mantle up to a few hundreds of kilometers 
depth. Density should increase in the top of the lithosphere, but may be followed by a den- 
sity drop where the second harzburgite layer is present. This is an unstable situation, but 
may ’freeze’ in (Nolet et al, 1986b). 
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Hypothesis testing 


Both the tectosphere and lithospheric doubling concept can be used for interpretation of the 
upper part of our models (z < 400 km). The tectosphere concept needs information about 
the depth at which there are no differences in mantle structure between different geological 
provinces, in our case the Scandinavian shield (SCSH1,2) and west-European platform 
(WEPL1,2). Lithospheric doubling looks for a strongly layered upper mantle and 
specifically an eclogitic layer that might mark the top of an obducted young ocean. Both 
models are not in conflict with each other like the pyrolite and piclogite models. 


5.3 Petrological Interpretation 


The assumptions in the hypothesis testing of the pyrolite and piclogite model, that acoustic 
velocities can be projected to greater depths along adiabats is not valid for the upper 200- 
400 km. Therefore Nolet et al. (1986b) parameterized a geotherm for this region and relied 
on laboratory measurements of seismic velocities at room temperature and pressure and 
their isothermal and isobaric derivatives with respect to pressure and temperature. By pro- 
jection along the adopted geotherm they were able to construct velocity curves for the most 
important mantle minerals. Although the laboratory measurements of the derivatives with 
respect to pressure and density are not very accurate, especially concerning the derivatives 
with respect to pressure and density, an important conclusion was reached: temperature and 
pressure effects cannot by itself cause high gradients in v, and p. As a consequence we must 
interprete large velocity or density gradients in terms of a change in mineralogy, chemical 
composition or partial melting. In figure 5.1.1. the upper 450 km of the WEPL1 and 
WEPL2 models are compared. 

The shear velocity contrast between the high velocity lid and the low velocity zone is 
enlarged in WEPL2 with respect to WEPL1. The high velocity lid in WEPL2 is confined to 
the depth range from 80-140 km. The corresponding density gradient indicates an increase 
in density, The average density in the 80-140 km depth interval amounts to 3.45 gr/ccm and 
can only be explained with the existence of garnet, since olivine, orthopyroxene and clino- 
pyroxene give density values of approximately 3.3 gr/ccm in this depth region (Nolet et al., 
1986b). This means that the interpretation of this high velocity/high density region as an 
eclogitic layer, as advocated by Nolet (1977) remains the most plausible one. No density 
inversion, as observed by Nolet (1977), is resolved at greater depth. 

The eclogitic layer near 100 km depth can only be explained by the lithospheric doubling 
concept: The basalt of the obducted ocean lithosphere is transformed to eclogite. A density 
inversion that marks the underlying harzburgite layer may not be resolved. 

The velocities in the low velocity zone are such that they do not coincide with any major 
mineralogy and partial melting is assumed. Since a possible Lehmann discontinuity is not 
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resolved, the gradient between 200 and 400 km depth is not very well determined. It can be 
observed however that both models, WEPL1 and WEPL2, exhibit the same gradient from 
250 to 350 km depth. 


density [gr/ccm] ; V, [km/s] 


= 1 | 1 ! 1 i A | 


“ 100.0 200.0 300.0 400.0 
depth [km] 
Figure 5.1.1, Density (lower curve) and shear velocity as a function of depth for 
WEPL1 and WEPL2. 


The deeper part of the model can be compared with shear velocity and density models for a 
pyrolite and piclogite composition by Bass and Anderson (BA ,1984) and Weidner (WE, 
1985). This is shown in figures 5.1.2 and 5.1.3. Concerning the shear velocity, it is obvious 
that the uppermost 350 km can not be fitted by any of these two models, which confirms our 
previous conclusions. At 400 km depth the jump in velocity agrees well with BA-piclogite 
and WE-pyrolite. A total disagreement exists with BA-pyrolite. The change in gradient 
between 500 and 600 km depth is compatible with WE-pyrolite, only at greater depth. This 
feature is not well resolved, but we are unable to maintain a gradient as in the BA- models 
in our inversion. Density is not well resolved below 150 km, but we can observe some 
features. The jump in density at 400 km depth is most likely to be larger than predicted by 
the BA-piclogite model. Inversions with a density difference smaller than WEPL2 generate 
a sharp density inversion after the 400 km discontinuity and this is thought to be unrealistic. 
At greater depths all models coincide with respect to density. 

Finally, we can compare our model for the scandinavian shield (SCSH2) with WEPL2., Both 
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models are displayed in figure 4.1.5, together with SCSH1. SCSHI1 and SCSH2 are very 
similar, except in the vicinity of a discontinuity and around 250 km depth where uncertainty 
exists related to the possible Lehmann discontinuity. 
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Figure 5.1.2. Shear velocity as a function of depth for WEPL2, Bass and Anderson 
(BA) piclogite and pyrolite and Weidner (WE) pyrolite. 


A general and well resolved difference between SCSH1,2 and WEPL2 is the low velocity 
layer at 160-220 km depth. This supports the tectosphere hypothesis that there are 
significant differences between shield and platform structure at least up to 250 km depth. 
The shield is much cooler, no partial melting seems to occur. In the upper part of the transi- 
tion zone shear velocity seems to be lower under the shield region, thus suggesting that 
there are differences in structure between shield and platform down to 500 km depth. This 
observation can also explain why PREM, as a average global model, exhibits systematically 
lower shear velocity values than our present inversion results. A difference in depth of the 
400 km discontinuity does not influence the shear velocity gradients in the transition zone, 
as shown in chapter 4. 

Regarding density we can only conclude that a possible density inversion exists under the 
Scandinavian shield between 200 and 350 km depth. This would also be in favour of the 
tectosphere model: the ’root’ is depleted in basaltic components and therefore less dense. 
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Figure 5.1.3. Density as a function of depth for WEPL2, Bass and Anderson (BA) 
piclogite and pyrolite and Weidner (WE) pyrolite. 


5.4 Conclusions 


The interpretation of shear velocity and density models, obtained by inversion of 
fundamental and higher mode Rayleigh wave phase velocities, has enabled us to investigate 
the structure and composition of the west-European platform. The most remarkable result is 
the existence of a layer of high velocity and high density, which is interpreted as an eclogite 
layer, at approximately 100 km depth. It is followed by a pronounced low velocity layer, 
which indicates partial melting. The transition zone is characterized by high shear veloci- 
ties and a change in shear velocity gradient at approximately 600 km depth. An interpreta- 
tion in terms of a pyrolite or piclogite structure is not possible, because of the uncertainties 
in elasticity parameters of the high pressure phases. However, the jump in density at the 
400 km discontinuity favours a pyrolite composition. The existence of an eclogite layer 
supports the lithospheric doubling theory. 

In the inversion procedure it became apparent that the existence of discontinuities in the 
starting model is of crucial importance for the inversion result. Introduction of these 
discontinuities, however, can lead to non-linear effects and as a result the linearized 
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inversion procedure is not valid any more. Therefore it is recommended to use the models 
produced in this study in a non-linear inversion procedure and test their validity as average 
models. A comparison with results from body wave modelling by Paulssen (pers. comm.) 
shows similarities, not only in the upper part of the model, but also regarding the high shear 
velocity values in the transition zone. 

The model for the Scandinavian shield is less reliable, because we compare a low frequency 
dataset for shield plus platform with a high frequency dataset for the platform only. As a 
result we can only look at gross differences. A major feature like the low velocity zone is 
clearly not present under the shield structure, while density seems to be lower under the 
shield from 200 to 350 km depth. In the transition layer there is evidence for a difference in 
shear velocity between shield and platform, suggesting a difference in structure between 
shield and platform down to 500 km depth. This observation is consistent with the tecto- 
sphere hypothesis. 
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APPENDIX A: 
The transfer function of a second order filter 


In this appendix the general form of a second order transfer function used in equation (1.3) 
is derived. 


Figure A.J. Circuit diagram of a second order Butterworth filter. 


For the voltages in figure A.1 we can write: 


i, in 
U;~Vo==>- U, -Vo= = 
Y, Y2 
- ig 
“ 0 ¥'; Y, ( ) 
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and for the current distribution: 


-ly=i;t13 


ing=ly 


Using (A.1) we can relate i3 to iz: 


and i, to i2 with the aid of (A.2) 
F 1 V3. 
pie yea 


Now we can use (A.1) to relate U,, to U; 


_YAtYoY3 
C2 0 = 
YiY2 
Substitution of i, 
. a S| 
lg=-l4= U 


finally leads to the response function: 


U, YY. 
U;  Yi¥AY,AY)+Y¥2t¥3) 


(A.2) 


(A.3) 


(A.4) 
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APPENDIX B: 
Optimization of a linear array using non-equidistant station location. 


1.0 


0.5 


O Kmin K max 


Figure B.1. Penalty function R(k) 


In stead of putting weights to array stations or developing more sophisticated estimators, a 
manipulation that can be carried out in the research center, one can vary station location 
provided the array is in its design phase. We will use a slightly modified version of the 
penalty function used by Cara (1978), now minimized with respect to station location: 


S(Aj,.....4,) =B cos6+ V sind 


Kein 
B= [RONG Pde; RO)=I ALP, kin 0k 
0 min 
Kaun 
V = [ Rk) IN()I? dk; R(k)=1, kin KmrinKmaxl (B.1) 
Kin 


where H (k ) stands for the array response function as a function of wavenumber k. Figure 
B.1 shows the penalty function R(k). If the interval of interest in higher mode analysis 
extends to k,,,,>k,, the aliasing wavenumber, and the influence of the mainlobe is limited 
to kin OME Can compute the penalty function 5. Minimization of this function S with 
respect to station location will provide a more optimal array response. 

For a better understanding of the behaviour of the penalty function we performed a simple 
excercise. In figure B.2 the array configuration has been modelled by means of two vari- 
ables x1 and x2, using x3 = f (x1,) to keep the total array length constant (In this case for 
NARS: 2600 km), 
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a simple exercise on choosing station locatim 


Configuration 


X)_.X2, Xs se Tepeated 
7 


i 2-3, 4 5.6 


[Ha] | i 


KEE SN 
SS SN 
ante 
Swe <—\ 
SN ES Se EN ON 
ASS SS 
SS SS 
SSSSES 


SS 
SESS 


Figure B.2. Penalty function S$(A,,...,A,) as a function of station 
Separation parameters x, and x2. For two minima and one maximum are the 
array response functions are plotted. 


Clearly a local maximum in S$ is observed for equidistant spacing, showing aliasing lobes 
as large as the main lobe itself within the interval of interest. Three local minima are visi- 
ble, each showing aliasing lobes removed from the specified interval. Furthermore we 
have been looking at a class of penalty functions in which 6 and k,,;,, were varied. The 
locations of the local minima didn’t change, although a trade off between the minimum 


value of S and k,;, exists. As 8 goes from 0 to = the minima become more pronounced. 


In conclusion we can say that the choise of non-equidistant station location in higher-mode 
surface wave analysis permits us to diminish the effects of spatial aliasing. 
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Appendix C: 
NARS event catalogue 


In this appendix we will present a list of all events recorded by the NARS array that are 
presently available on event tapes in GDSN format. In the last column the amount of 
stations that recorded these events is specified. These numbers include all triggers starting 
within one hour after the reported origin time. 


Time Latitude Longitude Depth mb nr of stations 


1983 0117 12:41: 38.092N 20.193E 
1983 01 24 08:17: 16.184N 95.149W 
1983 01 24 23:09: 12.952N 93.640E 
1983 01 26 16:02: 30.3558 179.427W 
1983 02 13 01:40: 39.945N 75.135E 
1983 02 14 03:20: 54.956N 159,.191W 
1983 02 26 07:10: 49.231N 155.617E 
1983 03 12 01:36: 36.5 4.103S 127.918E 
1983 03 18 09:05: 4.8578 53.507E 
1983 03 23 23:51: 06.3 38.333N 20.215E 
1983 04 03 02:50: 00.6 8.731N 83.115W 
1983 04 04 02:51: 34.8 5.730N 94.811E 
1983 0404 23:12:47.0 49.399N 155.640E 
1983 04 12 12:07: 54.5 4.8858 78.183 W 
1983 04 18 10:58: 47.8 27.721N 62.058E 
1983 04 30 14:03:49.6  41.523N 143.780E 
1983 05 01 18:10: 41.5 46.356N 153.490E 
1983 05 02 23:42: 37.8 36.242N 120.300W 
1983 05 02 23:42: 37.8 36.242N 120.300W 
1983 05 14 23:13: 46.3 38.488N 20.347E 
1983 05 15 00:24: 00.5 18.888S 175.710W 
1983 05 26 02:59: 58.8 40.480N 139.090E 
1983 06 09 12:49: 04.1 40.200N 139.050E 
1983 06 21 06:25: 27.5 41,315N 139.136E 
1983 0621 14:48: 5.5 24,.253N 122.480E 
1983 06 24 09:06: 46.6 24.197N 122.430E 
1983 06 24 07:18: 21.9 21.723N 103.380E 
1983 07 03 17:14: 23.1 9,.658N 83.644W 


& 


NAANDNWEHAUNMWWUAAKDATA TD WAWOAAIATDAAANAWA AN 
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Time 


nr of stations 


1983 07 05 
1983 07 11 
1983 07 12 
1983 08 06 
1983 08 17 
1983 08 25 
1983 09 07 
1983 09 12 
1983 09 21 
1983 10 04 
1983 10 17 
1983 10 22 
1983 10 28 
1983 10 30 
1983 11 16 
1983 11 24 
1983 11 30 
1983 12 02 
1983 12 22 
1983 12 30 
1984 01 01 
1984 02 01 
1984 02 01 
1984 02 07 
1984 02 10 
1984 02 11 
1984 02 16 
1984 03 05 
1984 03 06 
1984 03 19 
1984 03 24 
1984 03 27 
1984 04 06 
1984 04 20 
1984 04 22 
1984 04 23 


12:01: 27.3 
12:56: 28.0 
15:10: 03.3 
15:43: 52.5 
10:55: 52.8 
20:23: 32.5 
19:22: 04.8 
15:42: 08.2 
19:20: 42.8 
18:52: 12.8 
19:36: 21.3 
04:21: 36.3 
14:06: 06.3 
04:12: 27.6 
16:13: 00.0 
05:30: 34.8 
17:46: 00.4 
03:09: 05.6 
04:11: 29.2 
23:52: 40.5 
09:03:37.5 
07:28:27.8 
14:22: 07.6 
21:33: 20.5 
16:51: 21.0 
08:02: 50.0 
17:18: 42.5 
03:33: 51.1 
02:17: 21.0 
20:28: 39.8 
09:44: 02.5 
20:06: 33.3 
23:08: 23.2 
06:31: 10.1 
06:14: 21.6 
21:40: 35.5 


Latitude 


40.329N 
60.896S 
61.033N 
40.176N 
55.786N 
33.490N 
60.982N 
36.525N 
24.118N 
26.623S 
37.616N 
60.615S 
44.046N 
40.290N 
19.433N 
7.5498 
6.888S 
14.053N 
11.950N 
36.318N 
33.404N 
49.033N 
34.681N 
9.924S 
28.333N 
38.384N 
36.422N 
8.136N 
29.362N 
40.288N 
44.162N 
4.638S 
18.933S 
50.028N 
0.522N 
47 460N 


Longitude 
27,228E 
52.935W 
147.370W 
24.728E 
161.210E 
131.420E 
147.500W 
71.118E 
122.210E 
70.767W 
17.477W 
25.439W 
113.880W 
42.171E 
155.450W 
128.240E 
72.115E 
91.940W 
13.605W 
70.740E 
137.320E 
146.610E 
70.544E 
160.450E 
112.080W 
22.066E 
70.836E 
123.760E 
138.870E 
63.333E 
148.280E 
145.844E 
168.840E 
148.760E 
19.860W 
146,.720E 


Depth 


mb 


UANARUATRTTAZSOORAATUPWGZAWDAMNANDANADAHAANAAAG 
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Time 


1984 04 23 22:29: 
1984 04 24 04:11: 


1984 04 24 21:15 
1984 04 29 05:02 


58.3 
30.8 
: 19.0 
259.8 


1984 05 06 09:12: 01.6 
1984 05 07 17:49: 41.5 


1984 05 11 10:41 
1984 05 13 12:45 
1984 05 21 15:39 
1984 05 26 03:58 
1984 05 30 07:49 


1984 06 11 02:05: 
1984 06 15 14:22: 
1984 0621 10:43: 
1984 06 24 11:17: 
1984 06 24 14:30: 
1984 07 19 06:56: 
1984 07 19 23:25: 
1984 08 06 12:01: 
1984 08 06 19:06: 
1984 09 10 03:14: 
1984 09 14 23:48: 
1984 09 17 09:08: 
1984 09 18 17:02: 
1984 09 22 21:44: 
1984 09 28 00:03: 
1984 1015 10:21: 
1984 10 26 20:22: 
1984 1101 04:48: 
1984 11 20 08:15: 
1984 1202 06:09: 
1984 1203 04:08: 
1984 1228 10:37: 


49.8 
: 53.3 
: 00.5 
: 55.5 
243.5 
33.8 
23.0 
40.5 
12.0 
50.4 
11.0 
73 

54.1 
38.3 
08.8 
48.6 
52.6 
42.3 
18.3 
35.3 
07.5 
22.0 
49.8 
16.0 
44.0 
35.8 
53.6 


1985 02 28 20:53:47.8 


1985 03 01 17:11 
1985 03 02 15:47 
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117.5 
33.4 


Latitude 


22.044N 
30.792N 
37.320N 
43.274N 
38.817N 
41,766N 
41.832N 
42.927N 
32.671N 
43.8298 
4.8788 
30.722S 
15.7908 
35.349N 
17.993N 
36.853N 
52.912N 
28.131N 
0.1158 
32.396N 
40.346N 
35.761N 
32.0688 
33.973N 
32.012S 
25.795S 
15.7378 
39.176N 
8.147N 
5.214N 
20.423N 
44.237N 
56.174N 
27.462N 
2.082S 
1.964S 


Longitude 


99.147E 
138.360E 
121.690W 
12.571E 
25.661E 
13.886E 
13.948E 
17.731E 
121.500E 
39.142E 
151.600E 
71.214W 
174.860W 
23.291E 
69.345W 
3.768W 
4.198W 
129.530E 
122.520E 
131.800E 
126.850W 
137.480E 
178.370W 
141.440E 
178.390W 
176.010W 
173.780W 
71.340E 
38.758W 
125.240E 
115.750W 
148.140E 
163.570E 
128.449E 
119.670E 
119.727E 


175. 


129. 


202. 


mb 


5.9 
6.1 
5.8 
5.1 
5.0 
5.5 
5.3 
5.1 
5.6 
5.6 
6.1 
6.1 
6.1 
5.6 
6.0 
4.6 
48 
6.1 
6.1 
6.2 
6.1 
6.1 
5.6 
6.6 
5.6 
6.3 
6.5 
6.0 
6.5 
6.3 
6.0 
6.3 
6.1 
5.9 
5.7 
5.8 


nr of stations 
4 


TMMMAOANOAMWMWANGAWAAAWUAMOWAUAHAUMNEAWAMMN 


1985 03 03 
1985 03 04 
1985 03 04 
1985 03 05 
1985 03 09 
1985 03 13 
1985 03 16 
1985 03 17 
1985 03 18 
1985 03 19 
1985 03 25 
1985 03 28 
1985 04 03 
1985 04 09 
1985 04 19 
1985 04 19 
1985 04 20 
1985 04 23 
1985 04 24 
1985 04 30 
1985 05 01 
1985 05 02 
1985 05 06 
1985 05 09 
1985 05 10 
1985 05 14 
1985 05 14 
1985 05 15 
1985 05 16 
1985 05 19 
1985 05 19 
1985 05 20 
1985 05 31 
1985 06 03 
1985 06 03 
1985 06 06 


22:47:7.2 
00.32.21.8 
03:32:49.1 
13:40:10.2 
14:08:4.3 

19:34:57.6 
14:54:0.7 

10:41:38.4 
19:49:45.8 
04:01:8.0 

05:14:35.1 
16:07:06.8 
13:06:20.2 
01:56:59.4 
17:43:10.8 
17.55.34.6 
18:23:48.0 
16:15:12.0 
01:07:14.5 
18:14:12.7 
13:27:56.1 
08:55:16.3 
03:04:22.7 
19:05:21.5 
15:35:50.5 
13:24:57.8 
18.11.08.9 
20:12:45.8 
14:20:25.1 
08:07:48.2 
18.09.15.6 
15:11:40.6 
07:24:34.1 
02:45:32.0 
12.06.21.1 
02:40:12.9 


Latitude 


33.1358 
33.2078 
32.9258 
1.192N 
66.239N 
43.510N 
17.013N 
32.6338 
7.758N 
33.198S 
34.2548 
40.310N 
32.5848 
34.1318 
11.848N 
11.766N 
9.004N 
15.344N 
16.498N 
39.266N 
9.196S 
48.871N 
30.885N 
51.465N 
5.599S 
10.610S 
10.562S 
56.637S 
29.0818 
53.611N 
30.253S 
35.489N 
12.246N 
13.175N 
15,2898 
0.932N 


Longitude 


71.871W 
71.663W 
71.793W 
122.826E 
150.029W 
127,561W 
62.448W 
71.551W 
123.544E 
71.653W 
72.185W 
140.362E 
71.656W 
71.618W 
86.651 W 
86.851 W 
77. 460W 
120.610E 
120.815E 
22.810E 
71.230W 
156.329E 
70.269E 
177.913E 
151.045E 
41.423E 
41.424E 
25.330W 
77.735E 
160.526E 
71.329W 
87.173E 
144.280E 
90.138W 
173.516W 
28.432W 


Depth mb nr of stations 
33. 6.7 9 
33. 6.0 9 
33. 5.7 6 
33. 5.6 5 
12. 5.9 3 
10. 6.1 7 
13. 6.3 10 
33. 5.9 8 
33; 6.0 7 
42. 5.9 9 
45. 6.0 8 

166. 6.1 3 
33. 5.8 2 
38, 6.3 7 
72, 5.3 5 
60. 5.3 5 
38. 5.6 3 

188. 6.3 7 
33. 5.6 5 
27. 5.5 10 

600. 6.0 6 
43. 5.9 9 
37. 5.6 3 
33. 5.7 4 
27. 6.3 11 
10. 6.0 6 
10. 6.4 7 
33. 5.8 6 
10. 5.9 3 
63. 6.1 3 
39, 5.9 4 
33. 5.2 3 
32: 5.5 3 
66. 5.2 4 
33. 6.2 5 
10. 6.3 4 
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Time Latitude Longitude Depth mb nr of stations 
1985 06 10 15:37:03.3 28.1128 67.185W 180. 5.8 2 
1985 06 12 17:22:52.2 24.485N 122.078E 33. 5.2 3 
1985 07 03 04:36:52:4 4.463S 152.823E 41. 6.2 6 
1985 07 22 09:26:52.1 6.281S 148.725E 33, 5.9 5 
1985 07 28 22:59:54.5 60.3128 26.910W 33. 5.8 3 
1985 07 29 07:54:44.3 36.186N 70.890E 101. 6.7 3 
1985 08 02 07:46:51.4 36.173N 70.811E 103. 6.1 4 
1985 08 04 02:36:23.6 7.445N 123.463E 35. 5.8 4 
1985 08 11 00:19:02.4 11.157N 140.207E 33. 5.7 5 
1985 08 12 00:04:50.9 38.4208 73.490W 33. 5.5 5 
1985 08 12 03:49:17:9 37.739N 141.733E 51, 6.0 5 
1985 08 15 04:28:47.1 47.051N 18.063E 10. 48 5 
1985 08 21 11:26:28.8 9.2118 78.908W 61. 6.1 3 
1985 08 23 12:41:59.7 39.423N 75.274E 33. 6.4 9 
1985 08 28 20:50:48.4 21.0128 179.011W 625. 5.9 4 
1985 09 10 04:07:50.1 6.4508 149.888E 33. 5.7 5 
1985 09 11 20:45:51:7 39.328N 75.416E 33, 5.8 5 
1985 09 15 07:57:53.7 17.975N 97.183W 67, 5.9 2 
1985 09 19 13:17:47.8 18.182N 102.573W 33. 7.0 9 
1985 09 21 01:37:13.8 17.823N 101.671W 33. 6.3 11 
1985 09 26 07:27:48.9 34.6278 178.694W 33. 6.5 10 
1985 09 27 03:39:08.2 09.805S 159,.844E 30. 6.3 8 
1985 09 27 10:10:19:6 21.9448 174.804W 33. 5.8 3 
1985 09 27 16:39:48:6 34.507N 26.588E 59. 5.5 4 
1985 1005 15:24:02.2 62.257N 124.312W 10. 6.5 7 
1985 10 09 09:33:32.8 54.801N 159.573 W 31, 6.3 8 
1985 1012 22:20:37.6 0.857NS 29.876W 10. 5.3 8 
1985 10 13 15:59:53.5 40.317N 69.840E 33. 5.8 6 
1985 10 18 04:19:08.3 46.300N 146.287E 291. 6.0 2 
1985 1027 19:34:57.0 36.402N 6.746E 10. 5.5 8 
1985 1029 13:13:42.7 36.720N 54.805E 33, 6.0 4 
1985 1029 14:10:39:5 9564S 150.992E 10. 6.0 3 
1985 1031 19:33:07.1 53.258N 166.924W 33. 5.8 2 
1985 10 31 21:49:20:0 28.7478 63.186W 595. 5.8 3 
1985 11 07 19:12:29.8 35.1985 179.357W 33. 6.2 6 
1985 11 17 09:40:21.7 1.591S 135.012E 13, 5.9 6 
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1985 1121 21:57:14.9 
1985 11 28 02:25:42.5 
1985 11 28 03:49:54:1 
1985 12 14 06:46:13.3 
1985 12 16 02:44:35.6 
1985 12 16 08:04:06:1 
1985 12 21 01:13:21.0 
1985 12 23 05:16:03.9 
1986 01 15 20:17:31.4 
1986 02 03 20:47:35.8 
1986 0207 03:54:58.6 
1986 02 10 18:33:46.7 
1986 02 12 02:59:31.8 
1986 02 21 05:39:56.6 
1986 02 27 06:23:12.4 
1986 03 03 01:24:05.7 
1986 03 06 00:05:38.3 


41.720N 
13.976S 
13.9328 
3.651N 
11.719N 
14.153S 
14.0358 
62.207N 
21.376S 
27.889N 
45.408N 
39.521N 
36.352N 
43.346N 
24.069N 
41.945N 
40.393N 


Longitude 


19.320E 
166.209E 
166.167E 
126.699E 
85.854W 
166.356E 
166.508E 
124.302W 
170.284E 
139.402E 
27.865W 
143.186E 
140.932E 
25.954E 
122.303E 
20.265E 
51.534E 


nr of stations 


ay 


6 
6 
3 
6 
7 
8 
7 
8 
3 
3 
5 
3 
5 
3 
3 
5 
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APPENDIX D 
Dataset used in inversion. 


Rayleigh wave phase velocity measurements (in km/s) 


mode 0 
26/5/83 21/6/83 22/12/83 24/3/84 18/9/84 
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mean 


21/6+26/5/83__ 


mode 2 
T Is} | 21/6/83 22/12/83 6/3/84 24/3/84 ~—-18/9/84 ~—-21/6+26/5/83 mean _ st. dev. 
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106 


21/6/83 


22/12/83 


mode 4 


mode 5 


6/3/84 24/3/84 


18/9/84 mean st. dev. 


26.67 
23.70 
22.86 
21.33 
18.82 
17.78 
16.84 
16.00 
15.24 


6.44 
6.35 
631 
6.20 
6.03 
5.91 
5.82 
5.74 
5.70 
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Love wave phase velocity measurements (in km/s) 
for event 24/3/84 (84) and 18/9/84 (262). 


mode 0 mode | mode 2 mode 3 
84 262 84 262 84 262 84 | 262 
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Sn traveltime dataset used in inversion. 


Al[rad] | T[s 
.0524 82.4 


0698 
.0873 
1047 
1222 
1396 
1571 
1745 
.1920 
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130.4 
153.9 
177.6 
201.5 


225.5 
249.6 
272.8 
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